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ABSTRACT 

Limb-darkening is fundamental in determining transit lightcurve shapes, and is typ¬ 
ically modelled by a variety of laws that parametrize the intensity profile of the star 
that is being transited. Confronted with a transit lightcurve, some authors fix the pa¬ 
rameters of these laws, the so-called limb-darkening coefficients (LDCs), while others 
prefer to let them float in the lightcurve fitting procedure. Which of these is the best 
strategy, however, is still unclear, as well as how and by how much each of these can 
bias the retrieved transit parameters. In this work we attempt to clarify those points 
by first re-calculating these LDCs, comparing them to measured values from Kepler 
transit lightcurves using an algorithm that takes into account uncertainties in both 
the geometry of the transit and the parameters of the stellar host. We show there are 
significant departures from predicted model values, suggesting that our understanding 
of limb-darkening still needs to improve. Then, we show through simulations that if 
one uses the quadratic limb-darkening law to parametrize limb-darkening, fixing and 
fitting the LDCs can lead to significant biases -up to ~ 3% and ~ 1% in R p /R*, 
respectively-, which are important for several confirmed and candidate exoplanets. 
We conclude that, in this case, the best approach is to let the LDCs be free in the 
fitting procedure. Strategies to avoid biases in data from present and future missions 
involving high precision measurements of transit parameters are described. 
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1 INTRODUCTION 

It has been known since the first observations of our Sun’s 
surface that its observed intensity decreases towards the 
limb. This effect, termed limb darkening, crucially affects 
the shape of the transit signature a planet imprints in the 
observed stellar flux when passing in front of its host star. 
In practice, limb-darkening is parametrized by laws which 
depend on p = cos (#), where 9 is the angle between the line 
of sight and the normal to a given point of the stellar sur¬ 
face. Some of the most widely used limb darkening laws in 
exoplanet transit lightcurve fitting are given by 

= 1 — a(l — p) (the linear law), 

= 1 — ui(l — p) — M 2 (l — p) 2 (quadratic law), 
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= 1 — y' Cn (1 — /A 2 ) (non-linear law). 
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The non-linear law_ h as as a special case th e squ are-root 
law proposed by iDfaz-Cordovez fe Gimene3 (jl992l l , which 
is obtained by setting c 3 = C 4 = 0 , a nd the variant, three- 
parameter law, introduced by lSing et al.1 (12009I) . with ci = 0. 

The laws listed above do not include all 
the parametrizations availa ble for limb-darkening. 
iKlinglesmith fe Sobieskil (|l97Clh . for example, introduced 
the logarithmic law for early type stars, which is given by 


M 

10 ) 


— 1 — h(l — p) — hplnp, 


while [Cl are tl (l2003h introduced an exponential law given by 


= 1 - ei(l - p) - e 2 /(l - e M ), 


M. 
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but these laws are less used for transit fitting in prac¬ 
tice, mainly due to their complex forms which are harder 
to deal with computationally and are not implemented 
in the most widely u s ed transit modelling codes to date 
(|Mande l fc Agol 20021: Eastman. Gaudi fc AgoJ 201.il. see. 
howev er, Kiurkchieva et all ll201.lll : Abubekerov fc Gostevl 

jzoiah ). 
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Confronted with measurements of an exoplanetary tran¬ 
sit, observers usually model the effects of limb darkening 
by either fitting the limb-darkening coefficients (LDCs) of 
a given law or fixing some (or all) of them using tabu¬ 
lated values (see e.g., IClaret fc Bloemenl l201l| : ISind l20icil . 
for some recent tables using the Kepler bandpass). How¬ 
ever, exactly which strategy is optimal, an d in which situa¬ 
tions , is still unclear. Pre vious works (e.g., ICsizmadia et al.l 
21113 : I Muller et al ] sa have discussed this issue focusing 
on the uncertainty introduced on the parameters retrieved 
from transit observations but, to our knowledge, no study 
has yet addressed the potential bias introduced by them. 
Such study is called for as such biases could be limiting (1) 
the instruments currently obtaining high precision photom¬ 
etry like Kepler, (2) exoplanet population studies which by 
definition are based on averaging out ran dom errors but not 
systematic ones ('e.g.. ISchlaufmanl 120151 ) and/or (3) tech¬ 
niques that require high-precision measurements like trans¬ 
mission spectroscopy. Studying the biases introduced by the 
treatment of limb darkening and whether or not they are 
significant is the main aim of this work. 

The sources of potential biases introduced by assum¬ 
ing a given limb-darkening law can be separated in four. 
One issue is the differences in the tabulated values of the 
limb-darkening coefficients, which are large even when the 
same model ste l lar at mospheres are used. For example, 
ICsizmadia et al.l (120131 ') has shown that the differences be¬ 
tween different approaches at tabulating limb-darke ning co¬ 
efficie nts for the quadratic law with the ATLAS9 (IKuruca 
dzi) stellar atmosphere models can be as high as 20%, 
which if incorporated in the modelling can lead to signif¬ 
icant increases in the uncertainties of the retrieved tran¬ 
sit parameters ; _The second issue is the fact that, as shown 
bv lHowarthl (1201 lh . limb-darkening coefficients obtained di¬ 
rectly from the intensity profiles cannot be compared di¬ 
rectly to coefficients observed in transit photometry due to 
the fact that the optimization procedures are different in 
each setting. This in turn implies that using LDCs obtained 
from an intensity profile of the star as inputs in transit pho¬ 
tometry should lead to biases in the retrieved transit pa¬ 
rameters. The third issue is related to model complexity: 
we usually model limb-darkening, which in models is best 
represented by the non-linear law, with low-order laws such 
as the quadratic law. Finally, the fourth issue, and perhaps 
the most complicated to tackle, is the fact that we still do 
not know with certainty if available models do a good job 
at reproducing real intensity profiles of stars. 

The paper is organized as follows. In §2 we d etail our 
meth odology for obta ining LDCs from ATLAS9 (iKurucd 
GH) and PHOENIX (iHusser et al. ISoH) stellar model at¬ 
mospheres, and compare our results with the literature in 
order to try understand the discrepancy between previously 
published tables. In §3 we compare our model LDCs to Ke¬ 
pler estimates for a sample of star s, taking in consideration 
the effects mentioned by iHowarthl ([2011), and compare the 
performance of the ATLAS and PHOENIX models at pre¬ 
dicting the observed limb-darkening effect. In §4 we explore 
the biases introduced on different transit parameters by fix¬ 
ing or fitting the LDCs, §5 presents a discussion of our results 
and §6 the conclusions of our work. 


2 FITTING LIMB-DARKENING MODELS 

The fitting of the limb-darkening laws to intensity profiles 
obtained from stellar model atmospheres is, in principle, a 
relatively straightforward procedure. Given a normalized re¬ 
sponse function for a given telescope/detector S a (A) (e.g., 
f S\(X)dX = 1, although the normalization doesn’t really 
matter in practice for the calculation of limb-darkening co¬ 
efficients), one must integrate the specific intensity Ix{X, p) 
at each angle /i; = cos(#i) multiplied by the response func¬ 
tion, i.e., 

Hpi) = [ 2 Ix(X,pi)S x (X)dX, ( 1 ) 

J Ai 

where Ai and X 2 define the wavelength limits of the band. 
This gives the observed (by the instrument) intensity, which 
can then be fitted by any of the laws cited earlier after nor¬ 
malizing by 1(1). Of course, for CCDs, the recorded quantity 
are photons and therefore we have to divide the integrand 
in equation © by hc/X, where h is Planck’s constant and c 
is the speed of light. 

In this work we make use of two widely used model 
intensity libraries to calculate LDCs: the ATLAS9 model 
atmospheres, available from Robert L. Kuruc z’s webp agJi 
and the ID PHOENIX model atmospheres ([Husser et al. 
l2013f) . These models differ both in the geometry used to 
solve the stellar atmosphere (with ATLAS models using a 
plane-parallel approximation and the PHOENIX models us¬ 
ing a spherically symmetric atmosphere), and on the actual 
physics used on each of them (see, e.g., IPled l201ll . for a 
short up-to-date discussion on this matter), so a difference 
between the LDCs computed from both models is expected. 
We will derive and compare LDCs from these stellar at¬ 
mosphere models using the Kepler high-resolution response 
functior0 in order to compare our results with those in the 
literature. 

2.1 Fitting limb-darkening laws with the ATLAS 
models 

We first note that the ATLAS intensities are given per unit 
frequency and not per unit wavelength, so a a c/A 2 term 
has to be added in Equation ©, and so the integral in this 
equation now reads 

/(pi)= £ 2 SA(A)dA ’ 

where the limits of integration, Ai = 348 nm and A 2 = 970 
nm are, in our case, the first and last wavelengths present 
on the Kepler response function. We recall, however, that in 
order to fit the laws cited in the introduction what we want 
is the intensity normalized by /(l), i.e., 

J( M i) = Iv(X, pi)S\(X)/X dX 
f^K(X,l)S x (X)/XdX' 

1 http://kurucz.harvard.edu/grids.html 

2 http://keplergo.arc.nasa.gov/kepler_response_hires1.txt 
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Figure 1. ATLAS stellar intensity profile obtained for the Kepler bandpass for a G5V type star with solar metallicity (white points), 
with the colored lines showing fits to the most popular limb-darkening laws. In the upper panel, solid lines correspond to fits obtained 
following the method of S10, while the dashed lines correspond to fits obtained by following the method of CB11 (note these overlap in 
the leftmost panels). The lower panel shows the (percentual) difference between these fits. The vertical dashed black line marks // = 0.05 
(see text). 


After using numerical integration and interpolation to 
perform these integrals, the resulting normalized intensity 
profiles were fitted by a least-squares procedure using two 
different approaches in order to compare our results with 
previous works. The first approach was to fit all the angles 
for the non-linear law, while fitting only intensities at 
0.05 f or the rest of the laws, a procedure followed by ISind 
(boiol) . hereafter referred to as S10. The second approach 
was to fit to 100-points obtained by interpolating the IT 
angles given by the ATLAS modelfl in a linear grid from 
At = 0.01 to a i = 1.0 i n 0.01 steps using a cub ic spline, a 
procedure followed by IClaret fc Bloemenl (1201 lT) , hereafter 
refered as CBllQ. 

We note that, because the coefficients are linear with 
respect to the parameters in all of the laws considered, the 
least-squares solutions are unique and therefore have ana¬ 
lytic solutions given a set of intensities and angles. The exact 
solution, outlined in Appendix A, was used to perform the 
least-square fits. Figure [T] shows the intensity profiles for a 
G5V star of solar metallicity (i.e., with an effective temper¬ 
ature of T e g = 5500, log-gravity of log g = 4.5, [M/H] = 0.0 
and microturbulent velocity of r’turb = 2 km/s) where, for 
illustration, fits to the most popular limb-darkening laws 
used in the literature have been plotted following the meth¬ 
ods of S10 (solid lines) and CB11 (dashed lines). The lower 
panel shows the (percentual) difference between those two 
methods, illustrating that they actually differ by a small 

3 At first sight, fitting to 100 points that have been interpolated 
from a set of only 17 points given by the models seems to lack any 
justification. We show in what follows that there is a rationale for 
this procedure and that there is a well defined criterion under 
which this procedure gives better results. 

4 We will use S10 and CB11 to refer both to the papers and the 
methodologies assumed in those works to determine the LDCs 


(but, as we will show, significant) amount, with the differ¬ 
ences being larger for low order laws such as the linear (with 
a median difference of ~ 1%, reaching a ~ 2.5% difference 
near the limb) and smaller for high-order laws like the non¬ 
linear (with a median difference of ~ 0.1%, and reaching a 
~ 1% difference near the extreme limb). As one would ex¬ 
pect, the difference between the methods of S10 and CB11 
is more important for the low order laws. 

2.1.1 Comparing LDCs to previous results 

The large panels in Figure [5] show the limb-darkening co¬ 
efficients u\ and U 2 for the quadratic law, for a range of 
effective temperatures, TAj-, calculated by us using the AT¬ 
LAS models for stars with solar metallicity {[M/H] = 0.0), 
logg = 4.5 and = 2 km/s using the S10 and CB11 

procedures (blue and black points), along with the same 
calculations made by S1(Q (red points), and CB11 (green 
points). The small panels just below each large one depict 
the absolute difference between our values and the ones ob¬ 
tained by those previous studies. 

Overall, it can be seen that the coefficients we calculate 
follow similar trends as the ones presented in previous stud¬ 
ies, except for the region between T g g- « 7500 — 8500 (gray 
bands in Figure [5J where both our results and the ones of 
CB11 show a smooth decrease in the u i coefficient and a 
smooth increase in the u 2 coefficient, while the S10 results 
show a sharp increase and decrease in ui and M 2 respec¬ 
tively. We found that this is due to S10 using old versions 

5 The coefficients plotted here can 

be found in David Sing’s webpage: 

http://www.astro.ex.ac.uk/people/sing/David_Sing/Limb_Darkening.html 

In particular, we used the full version of Table 3. 
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Effective Temperature (K) 

Figure 2. Limb-darkening coefficients for the quadratic law ob¬ 
tained through the methods described in this work (blue and 
black), along with previous results by Sing (2010, red) and Claret 
& Bloemen (2011, green). Below each graph the absolute differ¬ 
ence between our results and each previous study is shown, while 
the difference between those studies is depicted by gray points 
and between the S10 and CB11 methods (but with our proce¬ 
dures) in black points. Note the large differences in the range 
T eff »s 7500 — 8500 (region denoted by the gray bands) with the 
Sing (2010) results; this is due to Sing (2010) using old versions 
of the ATLAS stellar model atmospheres. 


of the ATLAS stellar atmosphere^- Using those old models 
(available in Robert L. Kurucz’s webpage) we recover the 
overall shape of the coefficients presented by S10 but can¬ 
not, however, eliminate the differences between S10 and the 
LDCs we calculate using the S10 methods. 

Despite the similarity of the overall trend, there are 
significant differences between the LDCs we calculate and 
previous studies, even though we use their methods as de¬ 
scribed in their works and the same stellar atmospheres. The 
differences are of the same order of magnitude (~ 10%) for 
both coefficients of the quadratic law, with larger differences 
in the M 2 coefficient for cooler stars (T^g- < 5000 K), and in 
the Mi coefficient for hotter stars (T g g > 30000 K). Around 
solar temperatures (5000 K < T g g- < 6000 K) the differences 
are smaller, ~ 0.1 — 1% between our results and CB11, and 
~ 1 — 10% between our results and S10. Although we show 
here only the differences with the quadratic LDCs, there are 
very significant deviations between the coefficients of other 
limb-darkening laws too. In particular, we note that the 
limb-darkening coefficients of the non-linear law obtained 
by the different methods vary widely between works and, 
thus, these must be used with caution. We rule out that 
the differences arise from different interpolation and/or in¬ 
tegration methods by performing the same calculations in 
different ways and in different programming environments. 

Given the results above, it is clear that the differences 
with previous studies do not arise from different numerical 
approaches but, rather, they arise either from the actual 
method used to obtain the LDCs (i.e., from the very defini¬ 
tion of the integrals used to estimate the integrated intensi¬ 
ties) or from differences in the model atmospheres used to 
obtain those coefficients (e.g., older versions of the ATLAS 
model atmospheres). Unveiling the actual source of these 
discrepancies is currently not feasible as the actual codes 
used to obtain these tables are not publicly available. We 
believe that in the era of high precision measurements mak¬ 
ing all details available for scrutiny in order to try under¬ 
stand discrepancies is fundamental. Following this logic, we 
provide all the codes needed to reproduce the results in this 
paper through GitHutQ Our codes can be used to obtain 
LDCs for arbitrary response functions. 


2.1.2 Comparison of the methods used to fit the 
limb-darkening laws 

We now discuss which of the two fitting methods, i.e., that 
of S10 or that of CB11, is the most appropriate for obtaining 
LDCs using stellar model atmospheres, as they clearly show 
significant differences in Figure [2] (black points on the small 
panels). 

To assess which of the two methods is better, we define 
below a criterion that relies on an accurate analytic descrip¬ 
tion of the intensit y pr ofile given by the ATLAS models. To 
this end, we follow iHowarthl (l201ll ) and use the non-linear 


6 See the note on http://kurucz.harvard.edu/grids.html; the new 
models are the ones ending in *new.pck (e.g., im01k2new.pck). 
The old ones (e.g., im01k2.pckl9) are also available on the web¬ 
page (Kurucz, private comm.). 

' http: //™. github. com/nespinoza/limb-darkening 
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Figure 3. Differences between LDCs for the quadratic law ob¬ 
tained following the different methods discussed in this work (S10 
in red, CB11 in green) and the “real” underlying quadratic LDCs. 


limb-darkening law as an accurate descriptor of the inten¬ 
sity profile for the same stellar parameters used in Figure [2] 
We note that this was done not to try to exactly reproduce 
the model intensity profiles but, rather, to emulate an inten¬ 
sity profile that is similar to the one that would be actually 
observed in a real stellar atmosphere. From this profile, we 
then obtain quadratic LDCs by sampling N = 17 points at 
the same p values as the ones given by the ATLAS models 
and, with this, fit the intensity profiles using the methods 
of CB11 and a method similar to that of S10 (for this ex¬ 
periment, we choose to fit all the p values and not only the 
values for which p ^ 0.05, which is the original method of 
S10, in order to have a fair comparison between the coeffi¬ 
cients obtained by these two methods). 

In order to assess how good the retrieved coefficients 
are, we obtain the “real” underlying quadratic LDCs by 
defining them as the coefficients one would recover when 
sampling a very large number of points N from the (real) 
profile and then fitting those points with the quadratic limb- 
darkening law. Under this definition, when N —> oo, a set 
of “limiting coefficients” should be recovered, which are the 
“real”, underlying quadratic LDCs that best-fit the under¬ 
lying profile. We thus take the limit as N —» oo in our 
least-squares procedure, imposing that the data are sam¬ 
pled from an underlying profile described by the non-linear 
law (see Appendix B). The resulting limiting coefficients for 


the quadratic law are given in terms of the coefficients of 
the underlying non-linear law by 


w i,real 


lim u\ 

N —loo 


u 2 ,real 


lim M 2 

N oo 


12 164 

— Cl + C2 + T7^rC3 + 2c 4 , 
35 105 


10 


34 


-Cl — -C3 — C4. 

21 63 


(3) 


Figure [3] shows the results of this experiment. The fitting 
method of CB11 is clearly better at obtaining the “real” 
quadratic LDCs as we have defined them. The reason is 
that a cubic interpolation of the profile does a good job at 
retrieving the underlying profile (the non-linear law in the 
case of our experiment) and, thus, sampling points from this 
interpolation is very similar to sampling more points from 
the profile, and thus a result closer to the “real” coefficients 
follows. It is also interesting to see that both methods fail at 
retrieving the “real” mi coefficients for very hot stars. This is 
due to the fact that limb darkening for hotter stars is sharper 
than for cooler stars, changing abruptly towards the limb. 
Because the 17 angles given by the ATLAS models are more 
densely sampled close to the limb, this gives more weight to 
that region of the profile and thus the fit is dominated by 
it. Sampling uniformly across the profile, which is what the 
method of CB11 does, alleviates this problem and gives an 
overall better fit to the whole profile. 

As a final note on this topic, we would like to note 
that the limiting coefficients we have introduced are the best 
by construction because to obtain them we use a sampling 
scheme that weights the whole profile uniformly. To obtain 
the limiting coefficients one has to make use of the non-linear 
LDCs, which are obtained by fitting the model intensity pro¬ 
files using one of the discussed methods. However, as shown 
in the lower panels of Figure [lj in the case of the non-linear 
law the fitting method assumed is not very relevant, pro¬ 
ducing median offsets on the order of only 0.1%. Because of 
this, we believe that the “best” quadratic LDCs that one can 
choose are the “limiting coefficients”, whose input non-linear 
limb-darkening coefficients might be obtained from either of 
the S10 and CB11 methods, or some other similar sampling 
scheme. For practical applications, we recommend using the 
CB11 method to fit the non-linear law to a given intensity 
profile and then use the resulting coefficients to obtain the 
limiting coefficients according to equation ©• The codes we 
provide are able to carry out these calculations. 


2.2 Fitting limb-darkening laws to the PHOENIX 
models 

After performing the analysis for the ATLAS stellar atmo¬ 
spheres, we now describe how we obtain LDCs with the 
PHOENIX models. The method we used is very similar to 
the one described for the ATLAS models, except for the fact 
that now we have 78 angles and, unlike ATLAS models, the 
intensities are given per unit wavelength. There is an ad¬ 
ditional, very important difference: unlike previous studies, 
we are careful in taking in consideration that the geome¬ 
try of the PHOENIX models is different from that of the 
ATLAS models and, thus, the original intensity profiles ob¬ 
tained from these two model atmospheres are not directly 
comparable. 

Spherically symmetric models like the PHOENIX 
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ones extend the atmosphere by a small but impor¬ 
tant fraction (~ 0.4% for old er versions of PHOENIX, 
lAufdenberg. Ludwig fc Kervellall2004h from what one usu¬ 
ally calls the “stellar radius”, which is where the Rosse- 
land mean optical depth, tr, satisfies tr ~ 1. Plane-parallel 
model atmospheres like ATLAS, on the other hand, by def¬ 
inition have a very thin extended atmosphere outside the 
stellar radius and, therefore, have tr ~ 1 effectively at (or 
very close to) the limb (fx = 0). If we recall that the nor¬ 
malized stellar radius is given by r = y/l — fx 2 , this implies 
that for plane-parallel model atmospheres this outer “Rosse- 
land radius” is effectively at r = 1 (i.e., at /.i = 0), while 
for spherically symmetric models this radius is at r < 1 
(i.e., at /x > 0). This fact is evident when plotting the pro¬ 
files of plane-parallel and spherically symmetric model at¬ 
mospheres, where the PHOENIX models have a sudden de¬ 
crease in intensity around /x ~ 0.05 that is not present in 
the ATLAS profiles. This, of course, is not a fundamental 
difference between the models but, rather, a difference aris¬ 
ing from the fact that these two models are defining in a 
different way what the “outer stellar radius” or the “limb” 
is. 

One way of accounting for the different treatment of 
the stellar limb in spherical models is to search for the point 
t'max at which tr ~ 1 and re-define that point as having 
H = 0 before fitting any law to the intensity profile. This 
procedure allows to have meaningful and comparable pro¬ 
files between plane-parallel and spherically symmetric mod¬ 
els. According to empirical results with the PHOENIX mod¬ 
els, the point rmax in the intensity profile can be found by 
searching the value at which the derivative of the intensity 
profile with respect to the radi al intensity profile, \d l / dr|, 
is ma ximum (see, e.g.. IWittkowski. Aufdenberg fc Kervellal 
I2004I L Once found, it suffices then to re-normalize the radial 
profile by dividing by rmax, thus re-defining then the point 
at which r = 1 (p. = 0). 

Figured] shows the results of applying such corrections 
to the PHOENIX profiles for a G5V type star, where we 
have plotted the original PHOENIX profile, the “shifted” 
version of it and the ATLAS intensity profile for compari¬ 
son, which also illustrates the sudden decrease in intensity 
already mentioned for the PHOENIX models near p = 0.05 
(r ss 0.9987). The lower panel of Figure [I] shows, for illus¬ 
tration, the derivative of 7(r)/7(0), where the maximum is 
indicated by the dashed line. It is evident from the upper 
panel of Figured] that the (shifted) PHOENIX and ATLAS 
profiles agree for this star at the limb (which was expected 
because now both models are modelling the same portions 
of the stellar disk). It is also important to see the signif¬ 
icant changes between the original PHOENIX profile and 
the corrected ones: it is clear that the original PHOENIX 
models were not modelling the same portions of the stel¬ 
lar disk as the plane-parallel models and, thus, any LDCs 
derived from it couldn’t be directly compared to the AT- 
LAS models. We note that this correcti on was not made by 
IClaret, Hauschildt fc Wittel d2012l . 12013I ). hereafter referred 
as CHW, when deriving their “quasi-spherical” LDCs which 
were o btained in order to be com pared to ATLAS models, 
nor by iNeilson fc Lesteil d2013al fq). who also make direct 
comparisons between the spherical version of the ATLAS 
models and their plane-parallel counterparts. 

Figure [5] shows fits to the shifted PHOENIX intensity 


1 

• 1 '1 1 1 



— 0 — ATLAS model intensities 



—n— Original PHOENIX model intensities 

- 

0.75 

—□— Shifted PHOENIX model intensities 

- 

0.5 


- 




0.25 




1 







0 

_ 1 _1_ 1 _i_ 1 _1_ 

V. 



Figure 4. Close-up to the original (orange squares) and shifted 
(red squares) limb-darkening profiles for the PHOENIX model at¬ 
mosphere of a G5V type star with solar metallicity (top) and the 
derivative of the intensity profile as a function of r = Jl — fx 2 
(bottom). The black circles show the ATLAS profile for compari¬ 
son. The dashed line marks the value of rmax found in this case. 

profiles to the most popular limb-darkening laws used in the 
literature for the same model shown in Figure 2] following 
the methods of S10 (solid lines) and CB11 (dashed lines). 
Overall, it can be seen that in the case of the PHOENIX 
models a larger deviation is observed between the methods 
used to fit the profile than for the ATLAS models (compare 
the lower panels between this and Figure [TJ. In terms of 
following the profile, the CB11 method does a better job 
than the S10 method at following the whole profile in the 
case of the non-linear law. This was expected in light of 
our discussion in Section 2.1.2. As for the quadratic and 
linear law fits, both methods seem to do an adequate job at 
describing the profile given the low flexibility of the models. 

2.2.1 Comparing LDCs with previous results 

Figure [6] shows the analogue of Figure [2] for the quadratic 
LDCs obtained using the PHOENIX model atmospheres and 
the different methods discussed in this work, where a solar 
metallicity, logg = 4.5 and ^ ur b = 2 krn/s has been used. 
For comparison, the coefficients of CHW have been plotted 
in orange in the big panelfl We can see in the Figure that 
the overall trend for stars hotter than ~ 4000 K is similar 
for all cases except in the region between ~ 7500 — 9000 

8 We note that in the models that contain the intensity pro¬ 
files made available to the public by the PHOENIX team 
(http://phoenix.astro.physik.uni-goettingen.de/) there is 
no stellar model atmosphere for T e g = 5000 K, solar metallicity 
and log g = 4.5 and, thus, we have omitted this value published 
by CHW in this plot. 
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y = cos 8 


Figure 5. PHOENIX stellar intensity profile obtained for the Kepler bandpass for a G5V type star with solar metallicity (white points). 
The panels show the same information as Figure [Tj but for the “shifted” PHOENIX models. 


K, where the differences are on the order of ~ 50% and 
for stars cooler than ~ 4000 K, where the differences are 
larger than 100%. The overall differences in the coefficients 
are large in comparison to the ones seen for the ATLAS 
models, with median differences on the order of ~ 20% for 
both coefficients with previous works. This is mainly due 
to our approach using a new parametrization of the stellar 
disk for obtaining limb-darkening coefficients from spherical 
model atmospheres such as PHOENIX. 


3 MEASURING THE LIMB-DARKENING 

EFFECT FROM TRANSIT LIGHTCURVES 

A quantitative determination of the limb-darkening effect 
from transit photometry requires careful thought and has 
been t he subject of some scru tiny in recent years. For ex¬ 
ample, ICsizmadia et al.l d2013l) has shown that the presence 
of unnoculted spots can lead to significant changes in the 
LDCs, due to the fact that a spotted stars’ intensity distri¬ 
bution is more complex than the simple laws given in the in¬ 
troduction. Furthermore, even in th e case of unspot ted stars 
the interpretation of LDCs is subtle. IHowarthl (1201 il l showed 
that the geometry of the transit biases the LDCs obtained 
from photometry. In particular, he shows that high-impact 
parameter transits highly bias the observed LDCs because 
they sample chords of the stellar surface which are closer to 
the limb and, thus, sample a very different part of the inten¬ 
sity profile as the one sampled when fitting a whole model 
intensity profile. This is an unavoidable problem and, there¬ 
fore, LDCs obtained from high-impact parameter transits 
are not directly comparable to the ones obtained by model 
inte nsities. _ 

IHowarthl d201ll ) also stresses that, because the optimiza¬ 
tion processes are different when fitting an intensity profile 
directly from model stellar atmospheres than when fitting 
transit lightcurves, the coefficients obtained by those two 
procedures are not directly comparable even if the tran¬ 


sit is central. This is a very important and often over¬ 
looked fact: it implies that limb-darkening coefficients ob- 
tai ned from transit lig htcurves such as the ones obtained 
bv lMfiller et al.l d2013l) should not be compared directly to 
LDCs obtained from intensity profiles derived from model 
stellar atmospheres. That this is relevant can be verified in 
a straightforward fashion with a simple simulation study: 

(i) Select a good representation of a model intensity pro¬ 
file for a given set of stellar parameters derived from stellar 
model atmospheres, such as the non-linear limb-darkening 
law (i.e., select a set of coefficients ci, C 2 , C 3 , C 4 ). 

(ii) Generate a (noiseless) synthetic transit lightcurve by 
feeding the chosen representation of the model intensity pro¬ 
file of the star and using any set of geometric parameters for 
the transit. 

(iii) Fit this synthetic transit lightcurve with the same 
code as the one used to generate it, but now using a 
quadratic limb-darkening law parametrization, fixing all the 
geometric parameters of the transit to its input values (i.e., 
letting just the limb-darkening coefficients to float). 

The result of this experiment is always the same and, at 
first, counter-intuitive: the quadratic LDCs obtained from 
the input model intensity profiles (obtained through, e.g., 
the limiting coefficients obtained directly from the non¬ 
linear coefficients using equation ®, denoted in what fol¬ 
lows as (m, u 2 ), are always different from the ones obtained 
from the fit to the synthetic lightcurve, which we denote 
as { 111 , 112 ) ■ Figure [ 7 ] shows the results for this simple ex¬ 
periment, where we used a simulated transit lightcurve of 
a planet with a period of P = 3 days, semi-major axis to 
stellar radius ratio a/R „ = 0.1, planet-to-star radius ratio 
Rp/R * =0.1 on a circular (e = 0, uj = 0), edge-on (inclina¬ 
tion i = 7t/ 2) orbit (a typical “hot-Jupiter”). The non-linear 
law LDCs (ci, C 2 , C 3 and C 4 ) were obtained for stars with 
logg = 4.5, [M/H] = 0.0 and i^trirb = ^ km/s for differ¬ 
ent temperatures of interest for exoplanet studies using the 
ATLAS models via the method of CB11 and the Kepler re- 
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Figure 6. Limb-darkening coefficients for the quadratic law ob¬ 
tained through the methods described in this work (blue and 
black), along with previous results by Claret, Hauschildt & Wit- 
tie (2012,2013) in orange (here we plot the coefficients obtained 
with the “quasi-spherical models”, which only fit values of /z ^ 0.1 
using the original PHOENIX intensity profiles); below each graph 
the absolute difference between our results and this previous 
study is shown. 



U2 



Ui 


Figure 7. Simulation verifying the results of Howarth (2011), 
where the quadratic LDCs obtained by fitting the intensity pro¬ 
files from model stellar atmospheres (u \, U 2 ) for stars with log g = 
4.5, [M/H] = 0.0 and ^urb = 2 km/s using the ATLAS models 
for different temperatures are plotted against these same coeffi¬ 
cients (uJjMrj) obtained from synthetic transit light-curves gener¬ 
ated using these same intensity profiles. The gray line follows the 
temperatures of the simulated systems in increasing order for bet¬ 
ter visualization. The dashed lines depict the = u* lines, which 
the points should follow if the coefficients from model intensity 
profiles and transit light curves were directly comparable. 


sponse function, and from these we derived the limiting coef¬ 
ficients ( 111 , 112 )- The synthetic transit lightcurv e s wer e gen¬ 
erated using the formalism of iMandel fc Ago! (|2002l ), and 
the coefficients (u{, U 2 ) obtained via non-linear least-squares 
using the Levenberg-Marquardt algorithm. 1000 points be¬ 
tween phases —0.05 and 0.05 were generated in each simu¬ 
lated light curve. 

From the experiment we can see that, qualitatively, the 
effect of observing the LDCs from photometry is to overes¬ 
timate the “underlying” u\ coefficient, while the effect for 
the U 2 coefficient is to underestimate it. The exact map¬ 
ping ( 111 , 112 ) —> (u*,U 2 ), however, is non-trivial and, thus, 
the comparison between LDCs obtained from transit pho¬ 
tometry, ( ui , 112 ), and from intensity profiles obtained from 
model stellar atmospheres, (m, U 2 ), is also non-trivial. Fur¬ 
thermore, this mapping is in theory also dependent on the 
transit parameters, because the optimization process used 
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to fit a transit lightcurve is dependent on them, although 
we find this dependence to be usually negligible. There is 
thus no simple rule to compare LDCs obtained from transit 
photometry with the ones obtained by model intensity pro¬ 
files, and the comparison has to be done on a case-by-case 
basis. 


3.1 Comparing photometrically obtained 

limb-darkening coefficients to model values 

In order to compare the LDCs obtained from photomet¬ 
ric mea surements to va lues obtained from model intensity 
profiles, iHowarthl (|201ll i introduced the SPAM (Synthetic- 
photometry/atmosphere-model) algorithm, which can be 
summarised in three steps: 

(i) Fit a transit lightcurve and obtain the best-fit geo¬ 
metric parameters 9 P = {a/R*, R p /R,,e,io, i} and best-fit 
quadratic LDCs u* = {it{,it 2 }. Obtain the stellar param¬ 
eters that best represent the stellar host (e.g., via spectro¬ 
scopic observations) #» = {P e ff, log 5, [Af/fL], Ut ur b}- 

(ii) Generate a synthetic transit lightcurve using the best- 
fit geometrical parameters 6 P , and using an accurate repre¬ 
sentation of the model intensity profile of a star with stellar 
parameters #*, which we choose to be the non-linear law 
representation. 

(iii) Fit the synthetic transit lightcurve to obtain esti¬ 
mates of the geometric parameters and LDCs, but now us¬ 
ing a quadratic limb-darkening law to parametrize the limb- 
darkening effect. With this, obtain the coefficients u* = 
(uj, U 2 } which can now be directly compared to the best-fit 
coefficients obtained from photometry u* = {u{,ul}. 

There are three very interesting things to note about 
this algorithm. First of all, note that this algorithm actually 
tests the performance of the non-linear law, i.e., of one of the 
best representations of the model intensity profiles that we 
have available, because those coefficients are the ones that 
are used as inputs for the algorithm which are then com¬ 
pared to observations. This is very interesting, as it gives us 
a simple tool to asses how well our model atmospheres do 
at predicting real, observed, intensity profiles through tran¬ 
sit lightcurves. The second thing to note is that, because 
the non-linear law is the law being used as input, the ac¬ 
tual method for obtaining the LD Cs (i.e. , flux-conservation 
method, least-squares, etc., see iHowar thll201ll . for a discus¬ 
sion on the differences on the obtained coefficients by those 
methods for the case of the quadratic law) is irrelevant in 
practice, because all those methods give the same coeffi¬ 
cients in the case of the non-linear law, as shown bv IClaretl 
(1200011 . Finally, it is very important to note that this algo¬ 
rithm deals mainly with the bias introduced by modelling a 
complex intensity profile (e.g., one following a law close to 
the non-linear) with a simple, two-parameter law such as the 
quadratic one in the transit fitting procedure. This can be 
verified by noting that small changes in the input non-linear 
LDCs (e.g., obtained from stars with similar temperatures, 
metallicities or gravities) give rise to larger changes between 
the modelled (in) and photometrically observed (u*) LDCs 
than changes in the transit parameters. This means that 
although it is true that this algorithm deals with the bias 
associated with transit geometry, differences in this mapping 


between systems are currently mainly dominated by stellar 
parameters rather than by the actual transit geometry of a 
given system. 

Although a big step forward in how to properly compare 
limb-darkening coefficients, the SPAM algorithm has a few 
shortcomings: (1) the geometric parameters are fitted in step 
(iii), which is inconvenient as one would want to fix them, 
because the objective of the algorithm is to perform the 
mapping ( 111 , 112 ) (uf,uZ) given a geometric setting; and 
(2) the algorithm does not account for uncertainties on the 
geometric parameters of the transit lightcurve, 6 P , and/or 
for uncertainties on the stellar parameters, 9 *. The latter is 
a problem because, as stated above, changes in the geometry 
of the transit lead to changes in the limb-darkening coeffi¬ 
cients obtained from the SPAM algorithm, while errors on 
the stellar parameters can lead to significantly different coef¬ 
ficients ci, C 2 , C 3 , C 4 used in step (ii) to generate the synthetic 
lightcurves. Ideally one does not obtain only best-fit transit 
parameters for the geometric parameters, but also their pos¬ 
terior probability distribution functions (PDFs) given the 
data via, e.g., a Markov Chain Monte Carlo (MCMC) al¬ 
gorithm, which is information that we want to use. In most 
cases an MCMC approach for obtaining stellar parameters is 
not practical and therefore posterior distributions for these 
parameters are not usually available. One can still fit an ad- 
hoc distribution to a given set of stellar parameters and their 
respective errors (e.g., a Gaussian for the case of a parameter 
with symmetric error bars), and use that as an approxima¬ 
tion to the posterior distribution of the stellar parameters. 
We propose here a modified version of the SPAM algorithm 
that takes into account uncertainty information. We term 
this algorithm Monte-Carlo SPAM algorithm (MC-SPAM), 
and it consists of the following three steps: 

(i) Fit a transit lightcurve and obtain the posterior PDFs 
given the data for the geometric parameters, i.e., p(# p |data), 
with the corresponding posterior for the LDCs, p(u^|data). 
Obtain the posterior distribution of the stellar parameters 
(e.g., via spectroscopic observations), p( 0 »|data 2 LJ. 

(ii) Draw a set of geometric parameters 9 d and 
stellar parameters 9 d from the posterior distribution 
p(9 p , 0*|data, data 2 ) = p(0 p |data)p(#*|data2). Apply the 
original SPAM algorithm to these parameters and obtain 
u* ,d = {u*p d , the SPAM limb-darkening coefficients 
for the given draw. 

(iii) Repeat step (ii) to obtain a Monte-Carlo sample of 
the model/photometric LDCs u *, which can now be directly 
compared to the observed LDCs uL 

Because in most cases no MCMC chains are available 
for the stellar parameters, in order to sample directly from 
p(#*|data 2 ), we assume independence among the stellar pa¬ 
rameters and sample the parameters directly from their re¬ 
spective posterior marginal distributions, i.e., 


9 Note that we make explicit the fact that, in most cases, the 
dataset used to constrain the geometric parameters of the tran¬ 
sit and the dataset used to constrain the stellar parameters is 
different. The algorithm does not rely on this fact, however. 
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4 

p( 0 *|data 2 ) = JTp( 0 ;,*|data 2 ), 

i= 1 

where {6>i,*, 0 2 ,„ 0 3 ,*, # 4 ,*} = (T eff , log 3 , [M/H], 

For the modelling of the posterior marginal distributions, 
we assume the quoted estimates and errors come from an 
analysis done on a ~y surface, which is essentially an anal¬ 
ysis of the likelihood and, thus, can be approximated 

by an analysis of the posterior distribution of the param¬ 
eters, p(#*|data 2 ) ~ £(x 2 ), if one assumes flat priors for 
them. In practice, we model each marginal distribution by 
a Gaussian in the case of symmetric error bars, where the 
mean is set to the best-fit value of the parameter and its 
variance to the square of the error. For asy mmetrical erro r 
bars, we assume a skew-normal distribution llAzzalinilll985T) . 
which is the natural choice for a “normal-like” distribution 
with lack of symmetry. Details on the method we use to fit 
the parameters of this distribution given an estimate of a 
parameter and a set of asymmetrical error bars, as well as 
how to sample from the resulting distribution are given in 
Appendix C. We note that although log g, i^urb and T e g are 
positive quantities and our treatment can in principle allow 
negative values, in practice, the parameter uncertainties do 
not allow such values to be sampled and even if they did, 
our algorithm is capable of detecting and discarding those 
values as invalid. 

Our MC-SPAM algorithm is available at GitHulF’l. We 
now use it to compare our estimates of the observed limb- 
darkening coefficients to a set of Kepler planets. 


3.2 Comparing model to observed LDCs using 
Kepler data 

Photometrically derived LDCs from fits to Kepler transit 
lightcurves were obtained from various sources in the lit¬ 
erature in order to retrieve LDCs for stars with different 
parameters that host transiting planets with different ge¬ 
ometries. Because the posterior distributions for the param¬ 
eters of those systems are not published, we choose to use 
the same parametrization used for the stellar parameters in 
order to sample values given their published quantiles (i.e., 
in case of symmetrical errorbars we assume the posterior, 
marginalized distribution of each parameter is a Gaussian, 
while for the case of asymmetrical errobars we assume the 
posterior is best described by a skew-normal distribution). 
We note that the uncertainties associated with the transit 
parameters had negligible influence on the retrieved esti¬ 
mates of the limb-darkening coefficients obtained with our 
MC-SPAM algorithm, which were mainly dominated by the 
stellar parameter uncertainties. 

We took data from the high sig na l-to- noise, low-impact 
parameter sample of I Muller et all (120131 ) where we addi¬ 
tionally removed the objects that had impact parameters 
larger than b = 0.5, which were Kepler-43b ( b = 0.65), 
Kepler-45b (b = 0.6), Kepler-7b (b = 0.556), Kep ler-8b 
(6 = 0.72), KIC 5357901b fKOI-188b. iHebrard et aDl20l4 
b = 0.6), Kepler-41b (fo = 0.54) and Kepler-15b (6 = 0.56). 
In order to be as conservative as possible, we also removed 

10 http://ire#.github.com/nespinoza/mc-spam 


all the objects which have not been confirmed as planets 
to date, including Kepler-71b, for which no spectroscopic 
confirmation has been published so far that can rule out 
the po ssibility of it being a low-mass star, as iHowell et alj 
(|2010f ) can only constrain its mass to be less than 0.1 solar 
masses. We also remove Kepler-3b (HAT-P-llb) from our 
sample because its host sta r has been shown to have a sig¬ 
nificant amount of activity (iFraine e t al.ll20 14h and thus th e 
LDCs might be significantly biased ( Csizmadiajt_ah[ |20jjL 
In addition, we add Kepler-93b to our sample ( Ballard et all 
l2014ll . the recen tly validated Jupiter-si zed planets KOI-206b 
and KOI-680b (lAImenara et all 120151 1 and, at the expense 
of larger errors on the estimated limb-darkening coefficients 
but in order to expand the effective temperatures sampled 
in this work, we also add the newly confirmed planets with 
low impact-parameters Kepler-186f, Kepler-296f, Kepler- 
2966, Kepler-436b, Kepler-439 b, Kepler-440b, Ke pler-441b, 
Kepler-442b and Kepler-443b (I Torres et alJ l2015li . Table [T] 
summarizes the best-fit transit parameters for each of these 
planets, while Table [2] summarizes the host-star parameters. 
Table [3] presents both the model LDCs, (mi,M2), and the 
u* = (iq, U 2 ) coefficients, obtained using 1000 samples from 
our MC-SPAM algorithm for each target. We used both the 
ATLAS and PHOENIX models with the methods discussed 
in Section 2. A microturbulent velocity of 2 krn/s was as¬ 
sumed for stars for which no measurement of this parameter 
was published. 

Figure [8] shows the published limb-darkening coeffi¬ 
cients, (u{,« 2 ), as white points with errorbars. For eas¬ 
ier visualization, the planets have been divided into two 
groups: the ones with low precision (upper panels) and high 
precision (lower panels) limb-darkening measurements; note 
that these also separate the cooler (upper panels, T e g- = 
3572 — 5431 K) and the hotter (T e ff = 5520 — 7650 K) stars 
in the sample. At the sides of each datapoint, the MC-SPAM 
results, (mi,M;>) f° r each system using the ATLAS models 
(left, blue points) and PHOENIX models (right, red points) 
are shown. The arrows show the changes from the median 
of the model quadratic LDCs sampled by the MC-SPAM al¬ 
gorithm, (ui, U 2 ) (i.e., the values obtained directly from the 
intensity profiles; in this case, they are the “limiting coeffi¬ 
cients” obtained using the non-linear law), to the median of 
the resulting MC-SPAM values, (mi,U 2 ). Figure [9] shows the 
differences between the measured LDCs and the MC-SPAM 
values, with colored bands indicating the 68% band of the 
mean of those differences for the case of the ATLAS (blue 
bands) and PHOENIX (red bands) models. Kepler-296e was 
omitted from the calculation of the distribution of the mean 
of the differences (see below). 

In general, just as we observed in our simulation shown 
in Figure [7] the shifts between the (mi,M 2 ) and (m*,m 2 ) co¬ 
efficients produced by the MC-SPAM are larger for the 112 
coefficients than for the Mi coefficients, and the effect is to 
correct the model underestimation of the ui coefficients and 
the model overestimation of U 2 ; in other words, the effect of 
the mapping ui —» u\ is to increase the value of ui, while 
the effect of the mapping M 2 —> M 2 goes in the opposite di¬ 
rection. As observed in Figures [8] and [9] this mapping is in 
general very effective at better describing the observations, 
especially if the ATLAS model intensity profiles are used 
as inputs. It is interesting to note that the errors given by 
MC-SPAM are always dominated by the errors on the stel- 
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Figure 8. High (bottom panels) and low (upper panels) precision quadratic LDCs derived from transit photometry for several exoplanets 
(white datapoints) and MC-SPAM model LDCs (u^u^) using ATLAS (blue, to the left of each datapoint) and PHOENIX (red, to the 
right of each datapoint) models. The blue and red arrows next to the MC-SPAM results represent the mapping Ui —>• u*, i.e., from the 
original model LDCs obtained from fits to the intensity profiles (in practice obtained using the non-linear coefficients through equation 
and our MC-SPAM estimates. The temperature of the host star of each system is indicated above each of the planet names, inside 
the figures. Note the change in scale between the upper and lower panels. 
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Figure 9. Differences between the observed LDCs, (u^u^), shown in Figure [8] as white datapoints and the MC-SPAM estimates, 
using the ATLAS (blue) and PHOENIX (red) models, which are the blue and red points next to the white datapoints in 
Figure [8] The arrows represent the effect of the MC-SPAM algorithm in the observed difference between (u^u^) and the original model 
LDCs obtained from fits to the intensity profiles. The blue bands indicate the 68% bands around the mean of the differences using the 
ATLAS results (with the median of this indicated by the blue solid line), while the red bands show the same for the PHOENIX models 
(with the red solid line indicating the median of this distribution). These medians and the associated 68% values are also indicated next 
to each band (note that in the upper panels both bands overlap). The dashed black line marks zero, for reference. 
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lar parameters; this fact is especially evident for Kepler-77, 
whose uncertainty in the microturbulent velocity of ±0.3 
km/s leads to larger errorbars on the results using the AT¬ 
LAS models. This highlights the importance of estimating 
this parameter directly from spectroscopic measurements. 

From the sample of low precision LDCs, which is also 
the cooler star sample (upper panels of Figures [5] and O, 
one can see that, without taking into account Kepler-296e, 
the agreement with the data is very good for ATLAS and 
PHOENIX model atmospheres for the mi coefficients, which 
shows a mean difference of —0.046^if one uses the AT¬ 
LAS models and -0.049i°;i|* if one uses the PHOENIX 
models, both of which are consistent with a zero mean dif¬ 
ference. There is a barely significant offset of the 112 coeffi¬ 
cients (2.8cr for the ATLAS models, 2.9(7 for the PHOENIX 
models), with the model values apparently systematically 
overestimating those coefficients by a mean value of ~ 0.3. 
For the sample of high precision limb-darkening coefficients, 
which is also the hotter star sample (lower panels in those 
figures), the agreement of the ui coefficients is very good for 
the ATLAS models, with a mean difference of —0.008jlo'oi2 
which is consistent with zero, but poor for the PHOENIX 
models, whose mean difference is 0.083lQQQg, inconsistent 
with zero by more than 10cr. Note that this bias is very evi¬ 
dent and more prominent for stars hotter than 6000 K, where 
the model values overestimate the LDCs by ~ 0.1. For the 
U 2 coefficients, both model atmospheres show slightly signif¬ 
icant biases, with the PHOENIX models doing a better job 
at predicting the observed LDCs with a mean difference of 
—0.0241 q qjj , which is 2cr away from zero, and with the AT¬ 
LAS models showing a mean difference of 0.042+00*1, which 
is 3.8cr away from zero. 


As a final note, there is one particular object worth 
discussing in detail, Kepler-296e, which has LDCs which de¬ 
viate more than 2a from those of systems with host stars 
of similar st ellar paramet e rs, in cluding Kepler-296f, which 
according to I Torres et alj (120151 1 orbits the same host star 
in an orbit almost two times farther away from it. As we 
can see in Figures [8] and [9] the coefficient changes induced 
in Kepler-296e due to the geometry of the transit are not 
expected to be very different from those of Kepler-296f, so 
the geometry of the system cannot explain the differences 
on the observed LDCs. Activity could, in principle, produce 
significant biases on t he LDCs through unnoculted spots 
llCsizmadia et all 120130 . but it seems unlikely that activity 
affected only one of the observed transits. Becau se Kepler- 
296 is known to be a tight binary ijLissauer et al ■ 203 ) , one 
might be tempted to think that Kepler-296e mayb e did not 
orbit the same star as Kepler-296f as claimed bv lTorres et alj 
(l2015h . but its companion. However, both of the stars in the 
system have actually very similar spectral types and, thus, 
very similar LDCs, which implies that the observed LDCs 
are actually very different compared to both stars in the sys¬ 
tem. An analysis of other alternative hypotheses is out of the 
scope of this work, but we note that these are the kind of 
analyses that can be performed from measuring, comparing 
and interpreting LDCs from transit lightcurves using our 
MC-SPAM algorithm. 


4 THE EFFECT OF USING FIXED LDCS IN 

TRANSIT FITTING 

In the past se ction, we showed that, as first noted by 
iHowarthl (l201ll l. LDCs extracted from fits to intensity pro¬ 
files of stellar model atmospheres, (mi,M 2 ), are not directly 
comparable to the LDCs obtained from transit photome¬ 
try, This is due to the fact that the two optimiza¬ 

tion procedures are significantly different from one another 
and, thus, a geometry dependent mapping using synthetic 
lightcurves has to be carried out in order to obtain the coef¬ 
ficients (mi, M 2 ) that can then be compared to the observed 
LDCs. This implies that even if one could measure with ex¬ 
cellent precision the intensity profile of a given star, obtain 
its LDCs with that profile, and then measure those coeffi¬ 
cients from transit photometry, also with excellent precision, 
there is an expected bias between the two sets of coefficients. 
This in turn means that if one fixes the LDCs obtained from 
the intensity profile in the transit fitting procedure, then one 
is using potentially biased coefficients that can then lead to 
biased transit parameter^**!. A strategy to avoid this bias 
could be to let the LDCs as free parameters in the fit. How¬ 
ever, we also expect a bias on the transit parameters in this 
case if, as it is usually done, the intensity profile is modelled 
with a quadratic law which we we have seen is not able to 
accurately describe the full intensity profiles. 

I 11 addition to the above mentioned problems, there is 
the issue related to our imperfect knowledge of the under¬ 
lying, “true”, intensity profile. As we saw in §3, this is cur¬ 
rently an issue as our models are not able to reproduce the 
observed LDCs with sufficient accuracy. On top of this, ac¬ 
cording to our results in §2, there are differences even be¬ 
tween our own modelling of those profiles both between dif¬ 
ferent model atmospheres and between the different methods 
used to derive the LDCs from them. 

I 11 order to explore these sources of bias, in this section 
we perform simulations to study the possible biases intro¬ 
duced by our limb-darkening assumptions on the retrieved 
transit paramete rs, using transit l ightcu rves generated with 
the formalism of lMandel fc Ago] (120021 1. 

4.1 A simulation study 

In order explore the effect on the retrieved transit param¬ 
eters of fixing or having the the LDCs as free parame¬ 
ters, we simulate transit lightcurves with unit period, cir¬ 
cular orbits and an assumed intensity distribution for the 
host star of the transiting planet. The choice of units such 
that P = 1 is just for convenience of sampling directly 
in phase and has no consequences for what follows. The 
geometric parameters of the transit we can retrieve from 
our simulated light curves are the planet-to-star radius ra¬ 
tio, p = Rp/Rt, the semi-major axis to stellar radius ra¬ 
tio, or = a/R t , and the inclination of the orbit, i. The 
simulations were performed as follows. First, based on the 
data from all transiting planets discovered to date, we 
choose to generate synthetic transit lightcurves for plan¬ 
ets with all the combinations of parameters {a,R,p} with 

11 The level of bias will depend, among other factors, on the 
bandpass. In particular, issues noted in this paper should in gen¬ 
eral be less severe in the infra-red. 
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values a R = (3.27, 3.92, 4.87, 6.45, 9.52,18.18, 200} and p = 
{0.01, 0.06, 0.11, 0.16, 0.21}. In order to explore the effect of 
different impact parameters, b = cos (i)a R , we also varied 
b from 0 to 0.9 in steps of 0.1. This defines 350 different 
orbital configurations for our simulations. For each orbital 
configuration, we simulated 100 noiseless, uniformly sampled 
lightcurves with 1000 in-transit points and 400 out-of-transit 
points each, whose initial times were randomly perturbed. 
We first assume perfect knowledge of the underlying inten¬ 
sity profile by generating the transits using the non-linear 
law with the coefficients {ci, 02 , 03 , 04 } obtained in Section 2 
for models with logp = 4.5, solar metallicity, = 2 

km/s and effective temperatures between 3500 K and 9000 
K. We use the ATLAS models, the fitting method of CB11 
and the Kepler bandpass. Once we generate a simulated 
light curve, we retrieve its transit geometric parameters us¬ 
ing constrained non-linear least squares with the Levenberg- 
Marquardt algorithm using the lmfit packag43- We per¬ 
form the fit in two ways: ( 1 ) fixing the limb-darkening coef¬ 
ficients to the ones given by the quadratic law for the given 
star (using the limiting coefficients defined in Section 2); 
and (2) leaving the LDCs as free parameters in the fit. In 
the latter case, we fit for the parameters qi = (ui + U 2) 2 
and q 2 = Mi/ 2 (ui + uj), where the parameters q\ and q 2 
are constrained to be in ( 0 , 1 ), in order to obtain p hysically 
soun d solutions in our non-linear least-square fits (IKippingl 
l2013l h For each combination of the geometrical parameters, 
the median of the fitted parameters obtained from the 100 
generated lightcurves, along with the corresponding errors, 
are reported. 

4.2 The case of central transits 

Figure [10] shows the results of our simulations for the sim¬ 
plest case of a central transit, where the percentual bias in¬ 
duced on p and a R is shown as a function of the effective tem¬ 
perature of the host star used to extract the limb-darkening 
coefficients; the upper panel shows the bias induced when 
fixing the LDCs to model values and the lower panel shows 
the bias induced when letting those be free parameters in 
the optimization procedure (the retrieved inclinations in this 
case were not reported as they all arrived at the input val¬ 
ues). The sizes and colors of the points represent the input 
values of p and a R , respectively, with smaller, bluer points 
representing the lower values of both p and a R (p = 0 . 01 , 
a R = 3.27), and big, red points representing the highest in¬ 
put values of those parameters (p = 0.21, a R = 200). Note 
that the error bars are plotted, but are smaller than the 
smallest points in the plot. Note also that some points over¬ 
lap. 

We infer from our simulations that the biases are small 
for central transits (maximum of ~ 0 . 2 % on p in the case of 
fitting for the coefficients, ~ —0.05% when fixing them), 
although significant for several exoplanets: from a query 
done to the NASA Exoplanet Archival. 326 Kepler Ob¬ 
jects of Interest (KOIs) have quoted uncertainties on p lower 
than 0 . 2 % and for which the effects of this bias are impor¬ 
tant. Of those systems 57 are confirmed exoplanets and 113 

12 http://cars9.uchicago.edu/software/python/lmfit/ 

13 Query done on 2015/01/25. 


show quoted uncertainties lower than 0.05% (22 of which 
are confirmed exoplanets). It is interesting to note that the 
bias seems to be more important for deeper transits (big¬ 
ger points in the Figure), and the effect is to retrieve deeper 
transits and larger distances to the host star when fitting for 
the coefficients, while the opposite effect is introduced when 
fixing the coefficients, estimating slightly shallower transits 
and smaller distances to the host star than the real ones. 
Furthermore, the bias is clearly larger when fitting for the 
coefficients than when fixing them, which suggests that, if 
the underlying limb-darkening model is accurate at a higher 
level than the error made by fitting the coefficients, then the 
best strategy is fixing the LDCs to their model values. 

We note that the shape that the biases have as a func¬ 
tion of effective temperature are very similar between the 
different strategies employed to fit the transit lightcurves. 
For the Kepler bandpass, the difference between the Ui 
(LDCs obtained from model intensity profiles) and the u* 
(same coefficients recovered from transit fitting) observed 
in our experiment in the introduction of §3 (Figure [7| fol¬ 
low a similar shape with temperature as that observed for 
the biases, that is, starting at the cooler temperatures the 
offset is larger, decreasing until it reaches a minimum off¬ 
set at T e g- = 4750 K, then gradually increasing again for 
hotter host stars, slightly decreasing around T e g- = 8500 K. 
This means that the actual mapping m -4 u* is less severe 
for temperatures around T/g- = 4750 K, which is one of 
the reasons why we observe a smaller bias in the retrieved 
transit parameters around those temperatures in our experi¬ 
ments. Additionally, the performance of the quadratic law fit 
is also optimal around T g g = 4750 for this particular band¬ 
pass, becoming slightly worse for cooler and hotter stars, a 
fact already discussed in §2.1.2 (Figure0. The two points 
above serve to explain the observed shape of the biases with 
temperature and in particular the observed minimum of the 
biases at T e g = 4750 K. 

4.3 The case of low and high impact parameter 
transits 

Figures [lT] and [ 12 ] show the results of our simulations for the 
cases of low (b = 0.3) and high (6 = 0.8) impact parameters. 
We can see that the value of b strongly affects the observed 
bias in the retrieved transit parameters, showing biases as 
large as ~ 1% for p, and ~ 2% for a R and i. These biases 
are larger than for the case of central transits and, therefore, 
more important. A query to the NASA Exoplanet Archive 
returns 933 KOIs with uncertainties better than 1% on p, 
164 of which are confirmed exoplanets. Furthermore, the 
impact parameter not only modifies the order of magnitude 
of the observed bias, but also modifies the trends observed 
in the case of central transits. For example, in this case, the 
effect is more important for shallow transits than for deep 
transits at most temperatures. 

One very interesting fact about our simulations is that 
although for low impact parameters the bias seems to be 
larger for a R and i when fitting for the LDCs, for high im¬ 
pact parameters this effect is reversed, showing higher biases 
in those parameters when fixing the coefficients. The same 
effect can be seen for the bias in p in the case of small plan¬ 
ets around low tempera ture stella r host s. This is in agree¬ 
ment with the results of iHowarthl d201 il l. who showed that 
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Figure 10. Biases in the recovered planet-to-star radius ratio p = R p /R * and the semi-major axis-to-stellar radius ratio an = a/R „ as 
a function of temperature obtained from the simulations described in the text. The upper panels show results when fixing the LDCs, 
while the lower panels show the results when letting them to float in the transit light curve fit. The size of the points denote the input 
value of p (p = 0.01, small, p = 0.21 big points), while the color of the points represent the input value of an (an = 3.27, blue points, 
an = 200 red points). 


the difference between LDCs obtained from model atmo¬ 
spheres and the ones obtained from transit photometry in¬ 
creases as one increases the impact parameter of the transit. 
This implies that the fixed limb-darkening coefficient strat¬ 
egy should worsen as one increases the impact parameter, 
which is what we observe in our simulations. Therefore, for 
high impact parameter transits one should fit for the LDCs if 
one is interested in decreasi ng the bias, which contrasts with 
the suggestion of lMtiller et al.l |2013l ) of fixing them for high 
impact parameter transits based on the observed increased 
uncertainty in the retrieved transit parameters when fitting 
for the coefficients. In our view, the best strategy strongly 
depends on the quality of the photometry and the observed 
geometry: if the retrieved uncertainties when fixing the coef¬ 
ficients are smaller than the biases shown here, then the best 
strategy should be to let the coefficients be free parameters. 

4.4 The effect of an unknown stellar intensity 
profile 

In order to explore deviations from the known intensity pro¬ 
file, we make use of previously published limb-darkening co¬ 


efficients in order to simulate systematic offsets from the 
limb-darkening profile being modelled. To this end, we used 
the same transit lightcurves generated in the past sub¬ 
sections, generated with our non-linear limb-darkening co¬ 
efficients, but fitted them fixing the limb-darkening coef¬ 
ficients to the values published by CB11. In this way, we 
can explore the impact that different methods for obtain¬ 
ing the LDCs can have on the retrieved transit parameters 
which, we note, are only lower limits on the actual offsets 
between the modelled and real profiles, as both using our 
limb-darkening coefficients or the ones published by CB11 
in our analysis in §3.2 lead to offsets between observed and 
modelled LDCs of the same order. Figures fliil and fldl show 
the results of our experiments. 

As can be seen from our results, the biases follow a sim¬ 
ilar shape to the difference of the observed quadratic LDCs 
between our work and that of CB11 in Figure [5] that is, as 
expected, the biases are larger where the calculated LDCs 
show the largest differences. Also, the effect is again more 
pronounced for smaller values of p, with maximum biases on 
the order of 3% for planets around stars with T e gp = 4000 
K with moderate impact parameters (b = 0.3), temperature 
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Figure 11. Same as Figure ITOl but with an impact parameter of 0.3 and an additional panel showing the bias induced on inclination i. 
Note that the scale is different than that of Figure ITOl 
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Figure 12. Same as Figure [TT] but with an impact parameter of 0.8. Note that the scale is different than that of Figure [10] but the 
same as that of Figure [lT] 


and geometry at which we also observe biases on hr on the 
order of ~ 10%, as well as biases on the order of 5% on 
i. A query to the NASA Exoplanet Archive returns 2221 
KOIs with uncertainties better than 3% on p, 491 of which 
are confirmed exoplanets. This is a clear illustration of the 
importance of our limb-darkening assumptions. 


5 DISCUSSION 

5.1 Lessons learned from the comparison of 
model-to-observed LDCs 

In Section §3, we compared observed LDCs to model LDCs 
using our new MC-SPAM algorithm, which not only takes 
into account the geometry of the transit, but also the uncer¬ 
tainties in the transit and host star parameters. From these 
results, it is was apparent that overall the MC-SPAM re¬ 
sults using the ATLAS models do a better job at predicting 
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Figure 13. Same as Figure llOl but with an impact parameter of 0.3 and using LDCs obtained from a different intensity profile (obtained 
from the work of CB11) to the underlying one (generated using non-linear LDCs from this work). Note that the scale is different from 
that of previous figures. 
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Figure 14. Same as Figure ll3l but with an impact parameter of 0.8. Note that the scale between that figure and this one is also different. 


the observed u\ LDCs for hotter stars than the PHOENIX 
models, with the tables slightly turned in the case of the 
U 2 coefficients. Both models seem to do a good job at pre¬ 
dicting the Mi coefficients of cool stars, while slightly over¬ 
estimating their U 2 coefficients. These results suggest that 
the non-linear law in general seems unable to predict the 
limb-darkening effect with the accuracy necessary given the 
quality of available data and, because this is the best repre¬ 
sentation that we have for the ATLAS and PHOENIX model 
atmospheres, this suggests that those models are still not re¬ 
producing the observed limb-darkening effect with sufficient 
accuracy. As such, these model results have to be taken with 
care and one must be careful in trusting too much the LDCs 
derived directly from model atmospheres. Our view is that 
one should always be aware of the possible biases introduced 
by fixing those coefficients in practical applications. 

5.2 Limb-darkening and biases in transit 

parameters: implications for exoplanetary 
science 

The biases uncovered by our experiments in §4 are, in gen¬ 
eral, non-negligible. In particular, the large biases shown 
for the cases of non-central transits both when assuming a 
perfectly known intensity profile (§4.3) and when assuming 
deviations from it (§4.4) have several implications for exo¬ 
planetary science. First, the biases shown for p have a direct 
impact on population studies that rely on accurate measure¬ 
ments of the planetary radius an d their associated u ncertain- 
ties which, as recently shown by ISchlaufmanl J2015I ). impact 
the derived conclusions of these studies. It also has an impact 


on present and future studies that aim to constrain the in¬ 
terior composition of exoplanets from precise m easurements 
of the planetary radius and mass. For example, iDorn et alj 
J2015I) recently showed that a precision better than ~ 2 % 
on the planetary radius is one of the key ingredients in or¬ 
der to be able to constrain the interior structure of rocky 
exoplanets. In addition, the biases for an have an impact on 
derived quantities such as the calculated incident fluxes on 
exoplanets. This parameter, combined with the inclination 
i of the orbit, defines in turn the impact parameter, which 
according to our results can vary significantly depending on 
the different limb-darkening assumptions. 

5.3 Can the biases be corrected and/or avoided? 

As shown in in §4, the biases introduced by simple 
parametrizations of the limb-darkening effect such as the 
quadratic law are significant in comparison to the published 
parameters uncertainties for several confirmed and candi¬ 
date Kepler exoplanets. A natural question to ask is: how 
can one correct for those biases and how can one avoid them? 

The correction of biases in the transit parameters re¬ 
trieved from a given fit using the quadratic law would nec¬ 
essarily have to involve simulations such as the ones shown 
here around the parameters of interest. However, an even 
simpler method would be to fit the original lightcurve with 
higher order limb-darkening laws like the three-parameter or 
non-linear laws. The problem with this approach is that, al¬ 
though it seems obvious to switch to higher order laws once 
the order of magnitude of the bias surpasses the uncertain¬ 
ties in a given parameter, there is always an intrinsic bias 
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on the retrieved transit parameters related to our poten¬ 
tially imperfect knowledge of the exact shape of the, under¬ 
lying intensity profiles in stars which, as shown in §4.4, can 
lead to large biases if the coefficients are fixed in the tran¬ 
sit fitting procedure. As shown in Section 3, it appears that 
better modelling is indeed necessary given that the observed 
LDCs are not well reproduced by the models. This will be an 
important endeavour to minimize biases in the parameters 
estimated from transit light curves, and will be important 
to fully exploit data from both present high-precision mis¬ 
sions like Kepler, and future high-precision missions like the 
Transiting Exoplanet Survey Satellite (TESS) or the James 
Webb Space Telescope (JWST). 

Although a better modelling of stellar atmospheres is 
an important task, it is also a long-term one. Given the 
large amount of data already obtained by missions such as 
Kepler and the advent of future missions that will focus on 
retrieving high-precision transit parameters, it is important 
to discuss short-term strategies that can help overcome the 
biases due to our limb-darkening assumptions shown in this 
work. The most natural approach would be to fit the transit 
lightcurves with high order laws letting their parameters to 
float. However, this will most likely be a bad idea for low and 
medium signal to noise lightcurves due to the degeneracy 
between the parameters fitted in this procedure. One strat¬ 
egy to overcome this and the problems stated above could 
be to select a set of high-precision transit lightcurves with 
precisely known stellar host parameters that have enough 
signal-to-noise as to be able to be fitted with high order laws 
such as the three-parameter or non-linear law. This would 
in turn allow one to obtain fitted LDCs which could then be 
used to predict or constrain those coefficients for other stars. 
Although in theory those fitted LDCs are dependent on the 
transit geometry, this dependance as discussed in Section 3 
is rather weak in comparison with the typical uncertainties 
of the host star parameters and, thus, this strategy seems to 
be a promising one for achieving high accuracy and precision 
measurements of transit parameters. 


6 CONCLUSIONS 

In this paper, we calculated LDCs using the Kepler bandpass 
for both ATLAS and PHOENIX stellar atmosphere models, 
and showed that we cannot reproduce previously published 
values. We could not fully resolve the difference due to the 
fact the procedures used in the literature are not openly 
available, but we believe the differences between the differ¬ 
ent sources of limb darkening coefficients available in the 
literature are due to the use of different input model atmo¬ 
spheres. We make our codes available for users to reproduce 
our results and calculate limb darkening coefficients for any 
bandpass in a flexible way. 

We showed that different methods used in the fitting 
stage of the intensity profiles lead to different LDCs. We 
define a set of quadratic law limiting coefficients that are 
the best description of the “real” underlying intensity dis¬ 
tribution in the sense that they are obtained by sampling 
uniformly the whole profile. Among previously published 
methods using th e ATL AS models, the best is that of 
I Claret fe Bloemenl (120111 1 which by generating additional 
points to the ones directly available from the models via 


spline interpolation provides a good approximation to a 
denser sampling of the underlying profile. We also point 
out an important correction that needs to be applied when 
dealing with the spherically symmetric PHOENIX models so 
that the stellar radius definition is consistent with the plane- 
parallel ATLAS models. We provide updated limb darken¬ 
ing coefficients using PHOENIX models with this correction 
applied which are now directly comparable to those derived 
using ATLAS models. 

In order to map model LDCs to those determined by 
transit photometry we introduce an algorithm called MC- 
SPAM (Monte-Carlo Synthetic-Photometry/Atmosphere- 
M odel), whi c h bui lds upon the SPAM algorithm proposed 
by H owarthl (|201 il l . The algorithm takes into consideration 
not only the fact that LDCs obtained from intensity profiles 
of model stellar atmospheres are not directly comparable 
to the coefficients estimated from transit photometry, but 
also the fact that all the stellar parameters and geometrical 
parameters of the transit have measurement errors. We use 
MC-SPAM to compare limb darkening coefficients predicted 
by models for a sample of systems which have had their limb 
darkening coefficients determined from Kepler transit pho¬ 
tometry. We show that the ATLAS and PHOENIX models 
are not able to fully reproduce the observed limb-darkening 
effect for systems with a wide range of stellar parameters. Fi¬ 
nally, we showed that, when using the quadratic law, fixing 
and letting the LDCs be free in the transit fitting proce¬ 
dures induce biases on the retrieved transit parameters that 
are significant for several candidate and confirmed Kepler 
planets, even if one assumes a perfectly known stellar inten¬ 
sity profile. 

Given our results, we conclude that if one is confronted 
with a transit lightcurve and one decides to use the quadratic 
limb-darkening law to model the limb-darkening effect, the 
best strategy in order to minimize the bias introduced in 
the transit parameters is to let them float in the transit fit¬ 
ting procedure. Although a natural strategy to follow once 
the biases introduced by a given limb-darkening law are im¬ 
portant for a given system geometry and precision would 
be to switch to higher order laws, one has to be careful if 
the strategy is to fix the limb-darkening coefficients. This is 
because the strategy would also lead to an underlying bias 
due to our still incomplete knowledge of the intensity dis¬ 
tribution of real stars; our understanding of stellar model 
atmospheres is not good enough to avoid biases by switching 
to using fixed model coefficients from higher order laws. As 
a short-term solution to this problem, we propose to fit a 
set of high signal-to-noise lightcurves with precisely known 
stellar host parameters using a high order law parametriza- 
tion such as the three-parameter or non-linear laws in order 
to obtain the LDCs directly from data, which could then 
be used to predict or at least constrain the values of those 
coefficients for other exoplanetary systems. 

The results shown in this work imply that the current 
achievable precision is smaller than the current achievable 
accuracy if one considers only our limb-darkening assump¬ 
tions. Because of this, we call the attention to observers 
to be careful about their limb-darkening assumptions, and 
suggest always leaving room for flexibility on the LDCs in 
the transit fitting procedures even when fitting high order 
laws in order to shield against biases in the retrieved tran¬ 
sit parameters. If precision is more important than accuracy 
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(e.g., in applications such as transmission spectroscopy), we 
suggest performing simulations similar to the ones shown in 
this work in order to quantify the order of magnitude of any 
biases if the choice is to leave the limb-darkening coefficients 
fixed. 


7 ACKNOWLEDGMENTS 

N.E. is supported by CONICYT-PCHA/Doctorado Na- 
cional. N.E. & A.J. acknowledge support from the Ministry 
for the Economy, Development, and Tourisms Programa Ini- 
ciativa Cientffica Milenio through grant IC 120009, awarded 
to the Millennium Institute of Astrophysics (MAS). A.J. 
acknowledges support from FONDECYT project 1130857 
and from BASAL CATA PFB-06. We would like to thank 
the anonymous referee for helpful comments, questions and 
suggestions that helped to improve this work. We would also 
like to thank Pierre Kervella & Antoine Merand who pointed 
out to us the fact that PHOENIX models have a limb defini¬ 
tion inconsistent with that of plane-parallel ATLAS models, 
Benjamin Rackham for his assistance on writing initial ver¬ 
sions of the fitting procedures for the PHOENIX models, 
Antonio Claret, David Sing & Tom Evans for useful dis¬ 
cussions regarding the procedures they used for obtaining 
LDCs, and Ashley Villar, Jontahan Fraine, Amaury Triaud 
and Rafael Brahm for useful discussions regarding the re¬ 
sults of this paper. 

Some of the data presented in this paper were obtained 
from the Mikulski Archive for Space Telescopes (MAST). 
STScI is operated by the Association of Universities for 
Research in Astronomy, Inc., under NASA contract NAS5- 
26555. Support for MAST for non-HST data is provided by 
the NASA Office of Space Science via grant NNX13AC07G 
and by other grants and contracts. This paper includes data 
collected by the Kepler mission. Funding for the Kepler mis¬ 
sion is provided by the NASA Science Mission directorate. 


REFERENCES 

Abubekerov, M.K. & Gostev, N. Yu., 2013, MNRAS, 432, 
2216. 

Almenara, J.M, Damiani, C., Bouchy, F., Havel, M., Bruno, 
G., Hebrard, G., Dfaz, R.F., Deleuil, M., Barros, S.C.C., 
Boisse, I., Bonomo, A., Montagnier, G., Santerne, A., 
2015, A&A, 575, A71. 

Aufdenberg, J. P., Ludwig, H.-G. & Kervella, P., 2005, ApJ, 
633, 1, 424. 

Azzalini, A., 1985, Scand J Statist 12, 171. 

Ballard, S., Chaplin, W. J., Charbonneau, D., Desert, J.- 
M., Fressin, F., Zeng, L., Werner, M. W., Davies, G. R., 
Silva Aguirre, V., Basu, S., Christensen-Dalsgaard, J., 
Metcalfe, T. S., Stello, D., Bedding, T. R., Campante, 
T. L., Handberg, R., Karoff, C., Elsworth, Y., Gilliland, 
R. L., Hekker, S., Huber, D., Kawaler, S. D., Kjeldsen, H., 
Lund, M. N., Lundkvist, M., 2014, ApJ, 790, 12. 

Bonomo, A., Hebrard, G., Santerne, A., Santos, N., Deleuil, 
M., Almenara, J., Bouchy, F., Dfaz, R., Moutou, C., Van- 
huysse, M., 2012, A&A, 538, A96. 

Claret, A., 2000, A&A, 363, 1081. 

Claret, A. & Hauschildt, P. H., 2003, A&A, 412, 241. 


Claret, A. & Bloemen, S., 2011, A&A, 529, id. A75. 

Claret, A., Hauschildt, P.H. & Witte, S., 2012, A&A, 546, 
id. A14. 

Claret, A., Hauschildt, P.H. & Witte, S., 2013, A&A, 552, 
id.A16. 

Csizmadia, Sz., Pasternacki, Th., Dreyer, C., Cabrera, J., 
Erikson, A. & Rauer, H., 2013, A&A, 549, A9. 

Dfaz-Cordovez, J., & Gimenez, A., 1992, A&A, 259, 227. 

Dorn, C., Khan, A., Heng, K., Alibert, Y., Connolly, J., 
Benz, W. & Tackley, P., 2015, arXiv: 1502.03605 

Dunham, E. W., Borucki, W. J., Koch, D. G., Batalha, 
N. M., Buchhave, L. A., Brown, T. M., Caldwell, D. A., 
Cochran, W. D., Endl, M., Fischer, D., Furesz, G., Gau¬ 
tier, III, T. N., Geary, J. C., Gilliland, R. L., Gould, 
A., Howell, S. B., Jenkins, J. M., Kjeldsen, H., Latham, 
D. W., Lissauer, J. J., Marcy, G. W., Meibom, S., Monet, 
D. G., Rowe, J. F., Sasselov, D. D. 

Eastman, J., Gaudi, B. S. & Agol, E., PASP, 125, 83. 

Endl, M., Caldwell, D. A., Barclay, T., Huber, D., Isaac¬ 
son, H., Buchhave, L. A., Brugamyer, E., Robertson, P., 
Cochran, W. D., MacQueen, P. J., Havel, M., Lucas, P., 
Howell, S. B., Fischer, D., Quintana, E., Ciardi, D. R., 
2014, ApJ, 795, 151. 

Fortney, J. ,J., Demory, B.-O., Desert, J.-M., Rowe, J., 
Marcy, G. W., Isaacson, H., Buchhave, L. A., Ciardi, D., 
Gautier, T. N., Batalha, N. M., Caldwell, D. A., Bryson, 
S. T., Nutzman, P., Jenkins, J. M., Howard, A., Char¬ 
bonneau, D., Knutson, H. A., Howell, S. B., Everett, M., 
Fressin, F., Deming, D., Borucki, W. J., Brown, T. M., 
Ford, E. B., Gilliland, R. L., Latham, D. W., Miller, N., 
Seager, S., Fischer, D. A., Koch, D., Lissauer, J. J., Haas, 
M. R., Still, M., Lucas, P., Gillon, M., Christiansen, J. L. 
and Geary, J. C., ApJS, 197, 9. 

Fraine, J., Deming, D., Benneke, B., Knutson, H., Jorda, 
A., Espinoza, N., Madhusudhan, N., Wilkins, A., Todorov, 
K., 2014, Nature, 513, 526. 

Gandolfi, D., Parviainen, H., Fridlund, M., Hatzes, A. P., 
Deeg, H. J., Frasca, A., Lanza, A. F., Prada Moroni, P. G., 
Tognelli, E., McQuillan, A., Aigrain, S., Alonso, R., An- 
toci, V., Cabrera, J., Carone, L., Csizmadia, Sz., Djupvik, 
A. A., Guenther, E. W., Jessen-Hansen, J., Ofir, A., Telt- 
ing, J., 2013, A&A, 557, A74. 

Hebrard, G., Santerne, A., Montagnier, G., Bruno, G., 
Deleuil, M., Havel, M., Almenara, J.-M., Damiani, C., 
Barros, S. C. C., Bonomo, A. S., Bouchy, F., Daz, R. F., 
Moutou, C., 2014, A&A, 572, A93. 

Howarth, I. D., 2011, MRNAS, 418, 1165. 

Howell, S, Rowe, J., Sherry, W., von Braun, K., Ciardi, D., 
Bryson, S., Feldmeier, J., Horch, E. & van Belle, G., 2010, 
ApJ, 725, 1633. 

Husser, T.-O., Wende-von Berg, S., Dreizler, S., Homeier, 
D., Reiners, A., Barman, T. & Hauschildt, P. H., 2013, 
A&A, 553, id. A6. 

Kipping, D. M., 2013, MNRAS, 435, 2152. 

Kjurkchieva, D., Dimitrov, D., Vladev, A. & Yotov, V., 
2013, MNRAS, 431, 3654. 

Klinglesmith, D. A. & Sobieski, S., 1970, AJ, 75, 175. 

Kurucz, R. L., 1979, ApJS, 40, 1. 

Lissauer, J. J., Marcy, G. W., Bryson, S. T., Rowe, J. F., 
Jontof-Hutter, D., Agol, E., Borucki, W. J., Carter, J. A., 
Ford, E. B., Gilliland, R. L., Kolbl, R., Star, K. M., Stef¬ 
fen, J. H., Torres, G., ApJ, 784, 44. 


20 Espinoza & Jordan 


Lund, M., Lundkvist, M., Silva Aguirre, V., Houdek, G., 
Casagrande, L., Van Eylen, V., Campante, T., Karoff, C., 
Kjeldsen, H., Albrecht, S., Chaplin, W., Nielsen, M. B., 
Degroote, P., Davies, G., Handberg, R., 2014, A&A, 570, 
A54. 

Koch, D. G., Borucki, W. J., Rowe, J. F., Batalha, N. M., 
Brown, T. M., Caldwell, D. A., Caldwell, J., Cochran, 
W. D., DcVore, E., Dunham, E. W., Dupree, A. K., Gau¬ 
tier, III, T. N., Geary, J. C., Gilliland, R. L., Howell, S. B., 
Jenkins, J. M., Latham, D. W., Lissauer, J. J., Marcy, 
G. W., Morrison, D., Tarter, J., 2010, ApJL, 713L, 131 

Mandel, K & Agol, E., 2002, ApJ, 580, L171. 

Muller, H. M., Huber, K. F., Czesla, S., Wolter, U., 
Schmitt, J. H. M. M., 2013, A&A, 560, A112. 

Neilson, H. R. & Lester, J. B., 2013, A&A, 554, A98. 

Neilson, H. R. & Lester, J. B., 2013, A&A, 556, A86. 

Plez, B., 2011, J. Phys.: Conf. Ser., 328. 

Schlaufman, K., 2015, ApJL, 799, L26. 

Sing, D. K., Desert, J.-M., Lecavelier Des Etangs, A., 
Ballester, G. E., Vidal-Madjar, A., Parmentier, V., He- 
brard, G. & Henry, G. W., 2009, A&A, 505, 891. 

Sing D.K., 2010, A&A, 510, id.A21. 

Shporer, A., O’Rourke, J. G., Knutson, H., Szabo, G., 
Zhao, M., Burrows, A., Fortney, J., Agol, E., Cowan, N., 
Desert, J.M., Howard, A., Isaacson, H., Lewis, N., Show¬ 
man, A., Todorov, K., 2014, ApJ, 788, 92. 

Torres, G., Kipping, D., Fressin, F., Caldwell, D., Twicken, 
J., Ballard, S., Batalha, N., Bryson, S., Ciardi, D., Henze, 
C., Howell, S., Isaacson, H., Jenkins, J., Muirhead, P., 
Newton, E., Petigura, E., Barclay, T., Borucki, W., Crepp, 
J., Everett, M., Horch, E., Howard, A., Kolbl, R., Marcy, 
G., McCauliff, S., Quintana, E., 2015, ApJ, 800, 99. 

Wittkowski, M., Aufdenberg, J. P., Kervella, P., 2004, 
A&A, 413, 711. 



Table 1. Sample of confirmed Kepler planets for which limb-darkening coefficients have been obtained. 


Planet name 

Rp/R * 

i/° 

a/R t 


e 

id 

u{ 

u{ 

Reference 


Kepler-423b 

0 12662”*” 0 ' 00029 
u.izuuz_q .00028 

87.78l°-j° 

O 1 fi q+0.030 

o.lDo_o 034 

0 

(fixed) 

0 (fixed) 

0.4591°;°?? 

o.isil?;!?! 

Muller et al. 

12013) 

Kepler-77b 

0 09958”*” 0 ' 00026 
u.uyyoo_ 0 Q 0 Q26 

88.001°;“ 

9 749”*” 0,054 
»^»_0.055 

0 

(fixed) 

0 (fixed) 

0 514+ 0 - 019 
u - 01 -0.018 

0.1181?;??? 

Muller et al. 

12013) 

Kepler-17b 

0 13354”*” 0 ' 00017 

U.IOOd4_ 0 q 0019 

89.711°;?® 

5 707+ 9,9 ^ 9 
' u '-0.008 

0 

(fixed) 

0 (fixed) 

0 412”*” 0,015 

u.^iz_o 015 

0.1801°;!?? 

Miillcr et al. 

12013) 

Kepler-6b 

0 09412”*”°' 00008 
u.uy^i^_ 0< oooi3 

89.39l°; 2 ® 

7 551 +0.029 
' • ool _0.011 

0 

(fixed) 

0 (fixed) 

p 

Ip* 

oo 

o 

1 + 

O © 

o o 

o o 

-j -a 

0.1381°;!?? 

Muller et al. 

12013) 

Kepler-422b 

0 09625”*” 0,00025 
u - uyDZO _0.00022 

88.09l°;°6 

13.8261°;°?° 

0 

(fixed) 

0 (fixed) 

0-4381°;°! 

0 171 + 0,042 
1-0.047 

Miillcr et al. 

(2013) 

Kepler-12b 

fill C7Q+0-00013 

u.ii»/y_ 0 00013 

88.791°;“ 

o ni q+0.019 

o.Ulo_o 020 

0 

(fixed) 

0 (fixed) 

0-4281°;°°? 

0 124+ 0,018 
u - iZ4 _0.017 

Muller et al. 

12013) 

Kepler-2b 

0 07764”*”° ■ 00004 
U.U1 <O4_0 00004 

83.12l°;°s 

A 1 ro+0.006 

^±.ioz_0 006 

0 

(fixed) 

0 (fixed) 

H ‘357 - * -0 ' 007 

u.oo/ _ 0 007 

0 162"*” 0 ' 011 

U. 1DZ — 0.012 

Miillcr et al. 

12013) 

Kepler-5b 

n n , 7Q7o - l - 0-00006 
u.ui z,_ 0 Qooog 

89.39l°;°i 

6 459”*” 0 023 
u -^ jy _0.008 

0 

(fixed) 

0 (fixed) 

0.3681?;??? 

0 142+ 0 020 
0.019 

Miillcr et al. 

12013) 

Kepler-13b 

n ns^*^ 9,9 " 97 
u.uoooo_ 0 .00007 

85.821°;?° 

4.434l°;°?J 

0 

(fixed) 

0 (fixed) 

0.3081?;??? 

0 OOO+0-014 

u.^ -0.013 

Muller et al. 

12013) 

Kepler-93b 

0 014751 ”*” 0,000059 

U.Ul^/01_o 000059 

89.1831?;°?? 

12.4961°;°?! 

0 

(fixed) 

0 (fixed) 

0 44Q+ 0,063 
U. 2 ±‘4y_ 0 Q63 

H 1 oq+0.089 
U.loo_o 089 

Ballard et al, 

. (2014) 

Kepler-186f 

0 0205”*” 0,0012 
U.UZUO_o 0013 

89.96l°;“ 

1781!? 

0 

(fixed) 

0 (fixed) 

0 TO - * -0 ' 39 
u ‘ /u -0.38 

—0 1fi+ 0 - 31 
n-iO-o.30 

Torres et al. 

12015) 

Kepler-296f 

0 o362"*” 0 0022 
u.uooz_ 0 ooi8 

89.95l°;?2 

1371?? 

0 

(fixed) 

0 (fixed) 

0 35+ 0 - 31 
n.oO-Q 32 

—0 01+0.35 
U-U1 -0.20 

Torres et al. 

12015) 

Kepler-296e 

0 0297”*” 0,0029 
u.u^y/ _o 0037 

89.89l°;“ 

71 +29 
ri -10 

0 

(fixed) 

0 (fixed) 

1 e;o+0-32 
1 - OZ -0.34 

—0 70+ 0,23 
u - * u -0.27 

Torres et al. 

12015) 

Kepler-436b 

0 0354"*" 0,0024 
U.U9O^_0 0035 

89.93lg;S 

1041?? 

0 

(fixed) 

0 (fixed) 

0 67+ 0 - 44 

uo -0.48 

-0.141?;?? 

Torres et al. 

12015) 

Kepler-439b 

n no ,: iQ9+ 9 '"" 9 
u.uz,o»^_o 00111 

89.95l°;?2 

1421?? 

0 

(fixed) 

0 (fixed) 

0 39+ 0 - 28 
-0.34 

-0.081?;?? 

Torres et al. 

12015) 

Kepler-440b 

0.03038±°;?!?? 2 

89.93l°;°I 

99.0l?° 4 2 

0 

(fixed) 

0 (fixed) 

0.341?;! 

-0.051?;?? 

Torres et al. 

12015) 

Kepler-441b 

0.0280i°;°SiI 

89.971?;!? 

2601!? 

0 

(fixed) 

0 (fixed) 

0.281?;! 

—0 02”*”°' 30 
u - uz -0.20 

Torres et al. 

12015) 

Kepler-442b 

0.021li?;°°?® 

89.941°;°^ 

146l®° 

0 

(fixed) 

0 (fixed) 

1 f)2+ 0 - 50 
i.uo_o.49 

—0 43+0- 32 
n.^o_ 0< 43 

Torres et al. 

12015) 

Kepler-443b 

0 0304”*” 0 0022 
U.UoU4:_o 0022 

89.941?;°? 

15111 

0 

(fixed) 

0 (fixed) 

0 65”*” 0 ’ 41 

U.OO-048 

—0 1 5+°- 25 
u - iO -0.31 

Torres et al. 

12015) 

KOI-206b 

0.06590t°;°ool5 

89.2ll°;g 

6 44+0.62 
°‘^-0.62 

0 119+ 0 - 079 
U -- Lxy —0.079 

<»1m 

H qok+0.022 

u.oz,o_ 0 022 

0.2201?;??? 

Almenara et al. (2015) 

KOI-680b 

0.06384t°°°° 2 ° 

85-5il°;i 

6 35+ 0 - 51 
o.oo_ 0 51 

0 114”*” 0,077 
u ‘ 0.077 

1041?? 

0 374”*” 0 ' 024 
u.o/^_0 024 

0.1801?;!?? 

Almenara et al. (2015) 
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Table 2. Parameters of the host stars of each of the planets listed in Table [T| 


Star name 

W K 

log 9 (cgs) 

[M/H] a 

w turb/ km / s 

Reference 

Kepler-423 

57901^ 

4 57+ 0 - 12 
-0.12 

0.26l°;H 


Endl et al. (2014) 

Kepler-77 

5520t®° 

4 40+ 0 - 10 

0.20l°°| 

1.8 ±0.3 

Gandolfi et al. (2013) 

Kepler-17 

57811®* 

4 5^+ 0 - 12 
0.12 

0 26”*~°' 10 
U.ZO-Q 10 


Bonomo et al. (2012) 

Kepler-6 

5647 

4 236"*” 0,011 

0 34+ 0 - 04 

u -*5^_0.04 


Dunham et al. (2010) 

Kepler-422 

5972^ 

4 50 +0 ' 10 
^• ou -0.10 

rj 90+O.O9 
u.^o_o 09 


Endl et al. (2014) 

Kepler-12 

5947ti“ 

4 175+ 0 - 015 
0.011 

0-07l°;°l 


Fortne^et^aL (2011') 

Kepler-2 

6366 tH 

4 oi+°- 01 
^• u± -o.oi 

o-28i8;ll 


Lund et al. (2014) 

Kepler-5 

6297^gg 

3 96" 1 " 010 
o.»u_ 0 io 

0 04 + °- 06 
u,u4t —0.06 


Koch et al. (2010) 

Kepler-13 A 

7650llgg 

A 9+0.50 
^• z -0.50 

n 9+0.20 
U-Z —0.20 


ShDorer et al. (20141 

Kepler-93 

56691^ 

4 47Q+ 0 - 004 

4±.4t/u_ 0 004 

0 1 Q+0-10 

U.lO-010 


Ballard et al. (2014) 

Kepler-186 

37551®° 

4.7361°;°?^ 

—0 26”*~°' 12 
U.ZD_o 12 


Torres et al. (2015) 

Kepler-296 

35721*° 

4.833l°;°ll 

—0 12 - *” 0 ' 12 
u ‘ iZ -0.12 


Torres et al. (20151 

Kepler-436 

465i+ 100 

^ooi_l00 

4 619" 1 " 0,015 
^.ux»_o 028 

0 Qi+°- 10 

U-U± —0.10 


Torres et al. (20151 

Kepler-439 

543111°° 

4 514+0-035 
^• O14 -0.073 

0.02l°;l° 


Torres et al. (20151 

Kepler-440 

413411H 

4 706 - * -0 ' 049 

44. / u o—0.016 

-0.30l°:l® 


Torres et al. (20151 

Kepler-441 

43401111 

4 715~*~°' 047 

1 ±0 —0.024 

-o.57i8;ii 


Torres et al. (20151 

Kepler-442 

440211°° 

4 673+ 0 - 018 
^.u/o_ 0 021 

—0 37+ 0 - 10 

u '° -0.10 


Torres et al. (20151 

Kepler-443 

4723"*” 100 

4 614" 1 " 0,016 
^• ux ^-0.029 

—0 oi+°- 10 
U-U1 -0.10 


Torres et al. (20151 

KOI-206 

6360lH° 

q OQ9+0.056 

0.0^-0.056 

-o.oi±8:18 


Almenara et al. (2015) 

KOI-680 

6161111 

3 613" I_U ' U47 
o.uxo_o 070 

-o.i8i8:ll 


Almenara et al. (2015) 


a It is assumed that [M/H] = [Fe/H]. 
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Table 3. Results of the limb-darkening coefficients obtained using the MC-SPAM algorithm. 


Planet name 

U1 (ATLAS) 

u 2 (ATLAS) 

u\ (ATLAS) 

ul (ATLAS) 

U! (PHOENIX) 

u 2 (PHOENIX) 

u\ (PHOENIX) 

u ,I (PHOENIX) 

Kepler-423b 

0 41 q+0-023 
u.4iy_ 0 021 

0 263+ 0 ' 013 
u.zoo_o 014 

0 457 +O O 22 
u - 40 '- 0.020 

0 197+0-009 
u - iy ‘ - 0.010 

0 463+ 0021 
u. 400 — 0i011 

0 206+°'° 07 
u.zud_q 014 

0 504 - *" 0 021 
u.ou 4_ 0 013 

0 142+0-009 

u.i 4 z _ 0 014 

Kepler-77b 

0 367+ 0 ' 073 
u.oo / _ 0 12 8 

n 1 £9+0-039 

u.ioz-o 063 

0 394 + 0,077 
u.0U4_ o 138 

0 139+0-028 
u.ioy_ 0 047 

0 5Q3+0-0 ° 6 

u.uuo_o 005 

n 109 + 0.002 
U.XOZ-0 003 

0 546”*” 0 006 
u.O4o _ 0 006 

nil O+0.003 
U.XXO -0 003 

Kepler-17b 

0 420+°'° 18 
U.4tZU_o 016 

0 262+ 0 ' 010 
U.ZOZ_o oil 

0 463+ 0 ' 016 
u. 400 —Q.016 

0 190+°'° 06 
U.1»U_Q QQ 8 

0 464+ 0,017 
u. 4 O 4_ 0 009 

0 . 2061 °;°!° 

0 502”*" 0,017 

u.ouz _ 0 010 

0 145+0-007 

U.140_ooi2 

Kepler- 6 b 

0-448±°;°°® 

0 245+ 0 ' 006 
U.Z4to _ 0 005 

p 

rf*. 

00 

00 

1 + 
0 0 
b b 
0 0 

CO 00 

n 170 + 0.005 
U.l/(5_0 1 004 

0 495+ 0 006 

u.4tyo _ 0 0Q6 

n 1 Qc:+0-004 

U.XOO_Q QQ 3 

0.5381°;°°° 

OH 6+0-005 
U-ll°—0.004 

Kepler-422b 

0-384l°;°J6 

n 904 + 0.008 

U.zo4_o 010 

0 414+ 0 ' 016 
U.^i^_o.oi4 

n 995 +O OO 6 

U.^U_Q QQg 

0 446+ 0 ' 009 
u.^o _ 0 010 

n 91 c+0.005 

U.zxo_o 005 

0.4871!;!!! 

0 1 51 + 0.010 
U.xoi-o 007 

Kepler-12b 

n Q 70 + 0.016 
u.o < 0 —Q.016 

n 907 + 0.009 
u.z,o/ _o 009 

0 419+ 0 ' 014 
u .^±»_ 0 013 

0 206+°'° 04 
U.ZUO _ 0 005 

0 455+0-007 

u.4:00_o QOS 

0 209”*" 0 ' 003 
u -^uy _ 0 002 

0 496 -1 " 0,007 

U.^»U_Q QQg 

0 144+0- 003 
-0.003 

Kepler-2b 

0 324+ 0 ' 011 
u -°^_o.oio 

0 319+ 0 ' 005 
u.o ±»_ 0 005 

0 354+ 0 ' 012 
u. 004 —0.010 

n 940 + 0.003 
U.Z4o _ 0 003 

n 4 QQ+ 0.006 
U.4oo _ 0 005 

0 210”*” 0 ' 002 
u.zxu_o 002 

U.40D_Q QQ 3 

n 1 07 +O.OO 2 

u.xo/ _o 003 

Kepler-5b 

n 099 + 0.006 

u.oz,^,_o 006 

n qi 0 + 0.003 

u.o±o _ 0 003 

0 379 +O OO 5 

U.O 1 y — 0.005 

n 910 + 0.001 

u.z,±o_q ooi 

0 441 +o o ° 5 

U. 44 l_ Q 005 

0 209”*” 0 003 
u -^uy _ 0 001 

n 40 Q+ 0.005 
U.4oo _ 0 006 

o.i43i!;!!! 

Kepler-13b 

0 243+ 0 ' 045 
U.Z4to _ 0 021 

0 q^9+0-013 

u.ooz _ 0 028 

0 299+ 0 ' 038 
u.zyy _ 0 019 

0 246+ 0 ' 009 
U-Z4:U_q 012 

n QOI + 0.019 

u.ozi_q 012 

n 97 ^+ 0.010 

U.Z/O _ 0 Q 13 

n 971 + 0.023 
u.o/ x_o 013 

0.1981°;!!® 

Kepler-93b 

0 402+°'° 16 
u.mjz_o 017 

0 - 2681 °;°!° 

0 434 +O.OI 6 
-0.016 

n 91 9 +O.OO 8 

u.^x^_o 008 

0 465+ 0 ' 012 
U. 4 :OO _0 015 

0 20 C S ”*"°' 009 
u.zuo _0 007 

0 503+0 013 

u.ouo _0 016 

0 146+ 0 ' 011 
u.x 4 o _ 0 009 

Kepler-186f 

0-259l°;°37 

0-4201°;°^ 

n *31 q+0.039 

u.o±o_q 036 

0 303+°'° 31 

u.ouo_Q Q 28 

0 219+ 0 ' 035 
u.z,x»_o 035 

0 466 +0 ' 024 
u.^tou_o 033 

0 302+ 0 033 

u.ou ^_ 0 037 

o.309i!;!!! 

Kepler-296f 

0.289±g;g^ 

0 4Q4 4 " 0,025 
U.4U4 _ 0 02 g 

0 346+0-033 
U.04O_ o 032 

0 293+ 0 026 
u.zuo _ 0 018 

n 909 + 0.028 
u.zoz— 0 .027 

0 475+0 019 

u. 4 t/D_o 020 

n 099 + 0.027 

U.OZZ_Q Q 31 

0.3081°;!!! 

Kepler-296e 

0 293+ 0 036 

u.^»o _0 034 

0 401+°'° 27 
U -^ U1 _0.029 

0 345+0 0 32 
u.o^±o_o 036 

0 299+ 0 ' 031 

u.^»»_o 020 

0 234+ 0 ' 029 
u '^_ 0 , 02 9 

0 475+0 019 

u.4/0-0 022 

0 319+ 0 027 

U.ox »_ 0 Q 34 

0 316+ 0 ' 030 

u.oxo _ 0 016 

Kepler-436b 

0 617+0-013 

U.O1/-0 016 

011 5 + 0-011 
u.iio_ 0 > oi 0 

0 644+ 0 ' 015 

u.o44 _ 0 017 

n 071+0-Q14 
U-U< 1-0.012 

0 619+ 0 ' 011 

u.oiy _ 0 013 

nii 7+0.009 

U - 1J - 1 -0.007 

0 657+0-012 
u.od/_q 015 

0 057+0-01 2 
u.uu / _ 0 009 

Kepler-439b 

0 465+ 0 ' 022 
u.4too_o 021 

0 229+ 0 ' 014 
u.zzy _ 0 015 

0 493 + 0 022 
u.4tuo_o q 21 

0 179+0-014 

u.i/y_ 0 oi4 

0 505+ 9 ' 012 
u.uuu _ 0 012 

n 104 + 0.006 
0.1o4 _ 0 005 

0 545+0-012 
U.O4O _ 0 q 13 

0 1 iq+0.007 

U.liy_ 0> 007 

Kepler-440b 

0 41 q+0-107 
u.4iy _ 0 112 

n 977+0.092 
u.z/ 1 _q.088 

0.4721°;°°® 

O 1 77+0.068 
u>± 1 ‘ -0.062 

n 070 + 0.116 

U.o' o_o . 10 7 

n 9+0.083 

u.oxz-o 091 

0 441 +0-105 
u - 44i -0.092 

0 196+0-059 

u.iyo _ 0 063 

Kepler-441b 

0 479 + 0 ' 076 
»_ 0 120 

0 219+ 0 ' 090 
u.^,±»_o 058 

0.5201!;?!! 

0 146+ 0 ' 062 
u._mo_o 036 

0 456”*”°' 117 
u. z ±oo_ 0 12 g 

0 260 +0 ' 088 
u.zou_ 0 075 

0 507+°- 1 05 

U.OU/_Q 124 

o.i66i!;!!! 

Kepler-442b 

0 559 + 0 ' 032 
u.ooy _ 0 047 

0 159+ 0 ' 036 
u.iou_o 025 

0 590+°'° 29 

u.ouu _ 0 038 

0-1041“;°” 

0 546 +0 - 040 
u.o^o_Q 061 

O 1 77+0.046 
‘ -0.034 

0 5«5+°-037 

U.DoO_o 052 

o.iosi!;!!! 

Kepler-443b 

0 607+°'° 14 

u.ou^_ 0 018 

0 122+ 0 ' 014 
U. 1 ZZ -0 010 

0 633+0-016 

u.uoo _ 0 019 

0.0801°;°!! 

0 61 n+ 001 3 
u.°iu _ 0 014 

n 190 + 0.010 

U.XZO_Q QQ 8 

0 647+ 0 ' 015 

U.D4/-0 017 

0.0641°;!!! 

KOI-206b 

0 31 5+ 0015 

u.oio_o 013 

0 31 5+0-007 

u.o±o _ 0 006 

0 370+°'° 14 
u.o /u_o on 

n 91 c+0.005 

u.z,±u_o 003 

0 440+°'° 07 

U.44U_Q QQg 

0 209 +0 004 
u -^uy _ 0 oo 4 

O 4QQ+0-003 
U.400_Q QQg 

0 140+ 0 ' 004 
U.X4U_o 003 

KOI-680b 

n qoq+ 0 -Oio 
u.ozo_o 009 

0 3Q5+0-004 

U.oUO_o 004 

0 Q 5 Q+O.O 12 

u.ooo_o 010 

0 242+ 0 ' 004 
u.z^tz_o 006 

0 451+ 0 005 
u.^tox_o 003 

0 207” 1 " 0,001 

u.zu /_ 0 002 

0 503+0 005 

U.OUO_o 006 

n 1 90 +O.OO 6 

U. lzo _ 0 0Q4 
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APPENDIX A: LEAST-SQUARES FITS TO 
LIMB-DARKENING LAWS 

Fits of stellar model atmosphere intensities with the laws 
mentioned in the introduction require a least-squares proce¬ 
dure to be followed which, in general, minimizes the quantity 


N 

X 2 = ^2 Wi (y t - f(xi)) 2 , (Al) 

i=l 

where yi are the datapoints to be fitted, f(xi) the model 
for those datapoints, Wi the weight given to each squared 
residual (yi — f(xi)) 2 and N the number of datapoints to 
be fitted. In our case we set Wi = 1 and the objective is to 
minimize eq. dSB using the various models presented in the 
introduction of this work. 

If no constraint is provided, then the problem is easily 
solved for all the laws presented because the coefficients are 
linear in the parameters. Our objective then, is to minimize 
eq. ED with yi = I(fii)/I(l), and Xi = yi (the angles at 
which each of those integrations are done), for the differ¬ 
ent models f(yt) presented in the introduction, obtaining 
the optimal coefficients in a least-squares sense. The models 
have the form 

f([u) — 1 — 'y ] Ongn(yi), 


where 9 n are the parameters of the laws (e.g., 9 n = c„ 
for the non-linear law) and g n (yi) are the functions which 
make up each law, e.g., in the case of the non-linear law 
(which contains the linear, square-root and three-parameter 
laws depending on which coefficient 9 n one sets to zero), 
g n (yi) = (1 — yA 2 ). order to minimize eq. (IA1I) . we cal¬ 
culate the partial derivatives of \ 2 with respect to the dif¬ 
ferent coefficients 9 n and set them to zero. The calculation 
is easily found to give 

§£.X2(/M-7M//(1))^ = 0, 

i =1 


with 

df(yi) 

d9 k 


— 9k(yi), 


which, after rearranging terms, gives the system of k iinear 
equations for the n coefficients 9„ 


y ' J 9 n Oi n ,k — /3fc, k — 1,2, ..., 77., 

n 


with, 

N 

OLn,k — ^ ' 9n (^1)9 

i=1 
N 

Pk = 'y2g k (y i )(l- I(yi)/I( 1)), 

i= 1 

which are trivial to solve. Note also that if we write the linear 
system as A 9 = b, with A n}k = a n y, 9 = {#1, 9 2 ,..., 9 n } T 
and b = {/fi,/?2, ...,/3„} T , in this case the matrix A is sym¬ 
metric, so the system is not only linear but very fast to 
compute. As an example of how to use the above result, for 
the linear law the only function is gi(yi) = 1 — yi so, in this 


case, there is only one equation for the parameter 9\ = a 
with parameters 

N 

«1,1 = E( 1_ E 2 , 

i= 1 
N 

Pi = l]( 1 -M«)( 1 -/(Mi)A(l)), 

i= 1 


which gives, 



Ql.l 


Ef=r(i -^) 2 


(A2) 


APPENDIX B: LIMITING CASES FOR 
KNOWN TARGET LIMB-DARKENING LAWS 

If we try to fit an arbitrary law when knowing that we 
are sampling from a law of the form I(yi)/I( 1) (e.g., the 
non-linear law), then “limiting coefficients” corresponding 
to N —> 00 can be obtained. To perform this calculation, we 
sample N uniform yi points by defining yi = (i — 1)/( N — 1), 
with i = 1,2, ...,N (note that this samples yi angles from 
yi = 1 for i = N to yi = 0 for i = 1). We introduce this 
sampling sums that give the parameters (e.g., eq. (1 A2I) for 
the parameter of the linear law) and then take the limit as 
N —» 00. 


B1 Limiting coefficient a for the linear law when 
sampling from the non-linear law 


We sample N intensity (I(yi)/I( 1)) and angle (yi) pairs 
from the non-linear law, i.e., 

4 

I(yi)/I(l) = 1 - C4 1 - yl 12 ), 

n=l 


with known coefficients c n . We now fit the profile with a lin¬ 
ear law and determine the limiting coefficient a that follows 
as N —» 00. We first find the numerator and the denomi¬ 
nator in equation (IA2I) . The denominator is easily found to 
be 


~ ^) 2 = 


2 N 2 -N 
6 (N - 1) 


= h(N), 


while the numerator takes the form 

N 4 

E(i - ^X 1 - i)) = E 

i= 1 n=1 

with Sn = EZA 1 — k-i — Ei 12 + y ( y l+2)/2 ). It is straightfor¬ 
ward to show that 


S 1 


S 2 


S 3 


S4 


N 




1/2 


2 N 2 -N 
6(N — 1) 


h(N), 


M N N 

•W . 5/2 3/2 

y + -EE . 

i= 1 i= 1 

5 N 2 - AN 
12 (N - 1)’ 
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where we have expressed the sums of non-integer powers of 
Hi directly. Now equation (1 A21) reads 


E£i(1-A*)(1--TG*)/J(1)) 5, 

di — a = —■—*———:-;-= Cn 




h(N) 


and we need to take the limit as TV —> oo The main challenge 
for taking this limit is to obtain a closed form expression for 
the limit of the ratio between the sums of non-integer powers 
of Hi and h(N), terms which appear in the ratios Si/h(N) 
and Ss/h(N). To evaluate this, we first obtain an expression 
for the sums of non-integer powers of Hi- In order to do so, 
we note that, for (integer and non-integer) exponent k ^ ±1, 


E k 


l 

(TV — l) fc 
1 

(TV — l) fe 




. 3=0 




where in the last step we have used the Euler-Maclaurin 
summation formula. This implies that 


E*i^ 

h(N) 


6(7V — l) (1_fc l / N k 
(27V - 1) \k + 1 


^—+0(N k ~ 2 ) 


whose limit as TV —» oo is easily found as all terms to the 
right of the first term in this expression vanish in that limit, 
i.e., 


Etr/4 = lim 6(7V — i)( 1_fc ) 


lim 

n— nxj h(N) 


n— too (27V — 1) 

Using this result, we finally find 

4 

Sn - 

~ 10 


7V fc 

k + 1 


1 + k 


v-' 7 , 81 , 5 

lim a = lim > c„ = — ci + ci + — C3 + -C4. 

/ —loo AT—>no < ^ 


MTV) 


70 


B2 Limiting coefficients Mi and M 2 for the 
quadratic law when sampling from the 
non-linear law 


Following the results in Appendix A, in the general case 
of all the two-parameter limb-darkening laws there are two 
equations for the parameters 0i and 6*2, with 


0i 

02 


P 2 CX 2 .I — filOt2,2 
01 , 202,1 — 01 , 102,2 ' 
/7lOl,2 — /32Ql,l 
01 , 202,1 — 01 , 102,2 ’ 


(Bl) 

(B2) 


Specializing to the case of the quadratic law (i.e. 0i = u\ 
and 02 = U 2 ), we obtain 


N 

On.fc = £( l-^)" +fc . 

i=l 

N 

Pk = ^(l-M«) fc (l-7(/rO/7(l)), 

i=l 

for n = 1, 2 and /s = 1,2. We note that in this case a n ,k = 
afc.n, with 


Ol,l 


JV 


£u-mo 2 


27V 2 - TV 
6(7V — 1) ’ 


Ol 2 


= E('-^) 3 = 


AT 2 


i=l 

N 


o 2,2 = £(1 ~ m ) 4 — 

i=l 


4(7V — 1) ’ 

4 _ 7V(27V — 1)(37V 2 -37V- 1) 


30(TV - 1)3 


and the /3fc are of the form 


/3l — ^ ^ C n Sn-) 

n= 1 
4 

P 2 — ^ ^ CnB n , 


with the SVi given in the past sub-section and, 

7-2 »r JV 


Bi = 

B 2 = 
B 3 = 

B 4 = 


+ E2M l 3/2 -^ /2 -P l 1/2 , 


i=l 


27V — TV 
6(7V — 1) 

TV 2 

4(7V — 1)’ 

27V 2 — TV o 5 / 2 7 / 2 3/2 

6(7V - 1) + ^ ’ 

v y i=l 

97V 4 - 21A 3 + 147V 2 - TV 


30(7V - l) 3 


The denominator of the expressions for u\ and U 2 , i.e., the 
denominator of equations m and m, is given by 

2 (37V 2 — 37V + 2) (TV + 1)(2 — 7V)7V 2 

O = Oi 3 — Oi 1Q2 2 = 2 - 22 -nr 2 ---. 

i,2 12 2 ' 2 720(7V-1) 4 

We now note that the limits we want to obtain can be writ¬ 
ten as 

,. (7202,1 .. Pia. 2,2 

Inn Mi = hm-Inn -—, 

N— >00 N— >00 ot N—> oo Q; 

r /3iai,2 02Oti,i 

hm U 2 = hm-hm -—, 

N—> 00 N—>o o O'. N—> oo O 

which, using the same methods as in Appendix Bl gives 
12 164 

hm Ml = —Cl + C 2 + T-TT-C3 + 2c 4 
N-y oo 35 105 

10 34 

lm M2 = —Cl - —c 3 - c 4 

N—>0 O 21 63 


APPENDIX C: FITTING AND SAMPLING 
FROM SKEW-NORMAL DISTRIBUTIONS 
GIVEN PARAMETER ESTIMATES WITH 
ASYMMETRICAL ERRORBARS 

Given an estimate of a parameter 0 in the form 0E- 2 , where 
in general 04 ^ 02, we want to sample points from the poste¬ 
rior distribution of 0, given the data, only knowing that the 
distribution is asymmetric. One choice for performing such 
sampling is to assume that the distribution of 0 is a skew- 
normal distribution with parameters p, a and a, which is 
given by 

p{9\n,a,a) =p a 

where, 

p a (y) = 2<j>(y)$(ay), 
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and, 

4>{y) = exp(-y 2 /2)/v / 2^, $(ay) = f <j>(t)dt. 

J — OO 

In this distribution, y is the mean, cr 2 its variance and a, 
also called the shape parameter, defines its skewness. Note 
that when a = 0, we recover a normal distribution of mean 
H and variance cr 2 . 

Cl Fitting a skew-normal distribution to 
observed estimates 

A simple way to obtain the parameters {y, a, a} of the dis¬ 
tribution for each parameter is to assume that 6 — i T 2 , 6 
and 8 + ai define the 0.16, 0.5 and 0.84 quantiles of the pa¬ 
rameter distribution (i.e., the quoted value of the parameter 
defines the median and the errors define the 68% credibility 
bands around it). We can now easily formulate the prob¬ 
lem as a non-linear least-squares problem where the inde¬ 
pendent, observed, variables are the observed values at the 
given quantiles (i.e., x — {9 — 02 , 0, 9 + eri}) and the depen¬ 
dent variables are the quantiles (i.e., y = {0.16,0.5,0.84}). 
Then, we minimize 

r = || y-m(x,y,a,a) ||, 

where the i-th element of rh is given by 

/ Xi 

p(9\y,a, a)d9, 

-OO 

which can be solved with any non-linear least-squares algo¬ 
rithm such as Levenberg-Marquardt. 


C2 Sampling from a skew-normal with known 
parameters 

Once the parameters for a given skew-normal distribution 
are known, we want to sample values from it in a simple and 
efficient fashion. The following algorithm has been published 
by A. Azzalini in his personal webpag43, but we quote it 
here for completeness. 

First, one has to compute the parameter 8 = 

q/vT + a 2 . With this, one now samples the random vari¬ 
ables uo,v, both of which are independent and have stan¬ 
dard normal distributions, and generate the random variable 
Mi = Suo + %/l — 8 2 v which has correlation <5 with uo- Then, 
sample the random variable a which equals in if Mo > 0 and 
—Mi otherwise. With this, z has a skew-normal distribution 
with zero mean, cr = 1 and shape parameter a. Finally, the 
random variable SN = y + az has a skew-normal distribu¬ 
tion with mean y, variance cr 2 and shape parameter a. 

This paper has been typeset from a T^^K/ file prepared 

by the author. 
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